1 #ifndef lint 2 static char vcid[] = "$Id: posindep.c,v 1.5 1996/09/14 12:37:20 bsmith Exp bsmith $"; 3 #endif 4 /* 5 Code for Time Stepping with implicit backwards Euler. 6 */ 7 #include <math.h> 8 #include "src/ts/tsimpl.h" /*I "ts.h" I*/ 9 #include "pinclude/pviewer.h" 10 11 12 typedef struct { 13 Vec update; /* work vector where new solution is formed */ 14 Vec func; /* work vector where F(t[i],u[i]) is stored */ 15 Vec rhs; /* work vector for RHS; vec_sol/dt */ 16 17 /* information used for Pseudo-timestepping */ 18 19 int (*dt)(TS,double*,void*); /* compute next timestep, and related context */ 20 void *dtctx; 21 int (*verify)(TS,Vec,void*,double*,int*); /* verify previous timestep and related context */ 22 void *verifyctx; 23 24 double initial_fnorm,fnorm; /* original and current norm of F(u) */ 25 } TS_Pseudo; 26 27 /* ------------------------------------------------------------------------------*/ 28 /*@C 29 TSPseudoDefaultTimeStep - Default code to compute 30 pseudo timestepping. Use with TSPseudoSetTimeStep(). 31 32 Input Parameters: 33 . ts - the time step context 34 . dtctx - unused timestep context 35 36 Output Parameter: 37 . newdt - the time step to use for the next step 38 39 @*/ 40 int TSPseudoDefaultTimeStep(TS ts,double* newdt,void* dtctx) 41 { 42 TS_Pseudo *pseudo = (TS_Pseudo*) ts->data; 43 int ierr; 44 45 ierr = TSComputeRHSFunction(ts,ts->ptime,ts->vec_sol,pseudo->func);CHKERRQ(ierr); 46 ierr = VecNorm(pseudo->func,NORM_2,&pseudo->fnorm); CHKERRQ(ierr); 47 if (pseudo->initial_fnorm == 0.0) { 48 /* first time through so compute initial function norm */ 49 pseudo->initial_fnorm = pseudo->fnorm; 50 } 51 if (pseudo->fnorm == 0.0) *newdt = 1.e12*1.1*ts->time_step; 52 else *newdt = 1.1*ts->time_step*pseudo->initial_fnorm/pseudo->fnorm; 53 return 0; 54 } 55 56 /*@ 57 TSPseudoSetTimeStep - Sets the user routine to be 58 called at each pseudo-time-step to update the time-step. 59 60 Input Parameters: 61 . ts - time step context 62 . dt - function to compute timestep 63 . ctx - [optional] context required by function 64 65 @*/ 66 int TSPseudoSetTimeStep(TS ts,int (*dt)(TS,double*,void*),void* ctx) 67 { 68 TS_Pseudo *pseudo; 69 70 PetscValidHeaderSpecific(ts,TS_COOKIE); 71 if (ts->type != TS_PSEUDO) return 0; 72 73 pseudo = (TS_Pseudo*) ts->data; 74 pseudo->dt = dt; 75 pseudo->dtctx = ctx; 76 return 0; 77 } 78 79 /*@ 80 TSPseudoComputeTimeStep - Computes the next timestep for a currently running 81 pseudo-timestepping. 82 83 Input Parameter: 84 . ts - time step context 85 86 Output Parameter: 87 . dt - newly computed time-step 88 @*/ 89 int TSPseudoComputeTimeStep(TS ts,double *dt) 90 { 91 TS_Pseudo *pseudo = (TS_Pseudo*) ts->data; 92 int ierr; 93 94 PLogEventBegin(TS_PseudoComputeTimeStep,ts,0,0,0); 95 ierr = (*pseudo->dt)(ts,dt, pseudo->dtctx); CHKERRQ(ierr); 96 PLogEventEnd(TS_PseudoComputeTimeStep,ts,0,0,0); 97 return 0; 98 } 99 100 101 /* ------------------------------------------------------------------------------*/ 102 /*@C 103 TSPseudoDefaultVerifyTimeStep - Default code to verify last timestep. 104 105 Input Parameters: 106 . ts - the time step context 107 . dtctx - unused timestep context 108 109 Output Parameter: 110 . newdt - the time step to use for the next step 111 112 @*/ 113 int TSPseudoDefaultVerifyTimeStep(TS ts,Vec update,void* dtctx,double* newdt,int *flag) 114 { 115 *flag = 1; 116 return 0; 117 } 118 119 /*@ 120 TSPseudoSetVerifyTimeStep - Sets the user routine to verify quality of last time step. 121 122 Input Parameters: 123 . ts - time step context 124 . dt - function to verify 125 . ctx - [optional] context required by function 126 127 @*/ 128 int TSPseudoSetVerifyTimeStep(TS ts,int (*dt)(TS,Vec,void*,double*,int*),void* ctx) 129 { 130 TS_Pseudo *pseudo; 131 132 PetscValidHeaderSpecific(ts,TS_COOKIE); 133 if (ts->type != TS_PSEUDO) return 0; 134 135 pseudo = (TS_Pseudo*) ts->data; 136 pseudo->verify = dt; 137 pseudo->verifyctx = ctx; 138 return 0; 139 } 140 141 /*@ 142 TSPseudoVerifyTimeStep - Verifies that the last time step was ok. 143 144 Input Parameter: 145 . ts - time step context 146 . update - latest solution 147 148 Output Parameter: 149 . dt - newly computed time-step (if it had to shrink) 150 . flag - indicates if current timestep was ok 151 152 @*/ 153 int TSPseudoVerifyTimeStep(TS ts,Vec update,double *dt,int *flag) 154 { 155 TS_Pseudo *pseudo = (TS_Pseudo*) ts->data; 156 int ierr; 157 158 if (!pseudo->verify) {*flag = 1; return 0;} 159 160 ierr = (*pseudo->verify)(ts,update,pseudo->verifyctx,dt,flag ); CHKERRQ(ierr); 161 162 return 0; 163 } 164 165 /* --------------------------------------------------------------------------------*/ 166 167 static int TSStep_Pseudo(TS ts,int *steps,double *time) 168 { 169 Vec sol = ts->vec_sol; 170 int ierr,i,max_steps = ts->max_steps,its,ok; 171 TS_Pseudo *pseudo = (TS_Pseudo*) ts->data; 172 double current_time_step; 173 174 *steps = -ts->steps; 175 176 ierr = VecCopy(sol,pseudo->update); CHKERRQ(ierr); 177 for ( i=0; i<max_steps && ts->ptime < ts->max_time; i++ ) { 178 ierr = TSPseudoComputeTimeStep(ts,&ts->time_step); CHKERRQ(ierr); 179 current_time_step = ts->time_step; 180 while (1) { 181 ts->ptime += current_time_step; 182 ierr = SNESSolve(ts->snes,pseudo->update,&its); CHKERRQ(ierr); 183 ierr = TSPseudoVerifyTimeStep(ts,pseudo->update,&ts->time_step,&ok); CHKERRQ(ierr); 184 if (ok) break; 185 ts->ptime -= current_time_step; 186 current_time_step = ts->time_step; 187 } 188 ierr = VecCopy(pseudo->update,sol); CHKERRQ(ierr); 189 ts->steps++; 190 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol);CHKERRQ(ierr); 191 } 192 193 *steps += ts->steps; 194 *time = ts->ptime; 195 return 0; 196 } 197 198 /*------------------------------------------------------------*/ 199 static int TSDestroy_Pseudo(PetscObject obj ) 200 { 201 TS ts = (TS) obj; 202 TS_Pseudo *pseudo = (TS_Pseudo*) ts->data; 203 int ierr; 204 205 ierr = VecDestroy(pseudo->update); CHKERRQ(ierr); 206 if (pseudo->func) {ierr = VecDestroy(pseudo->func);CHKERRQ(ierr);} 207 if (pseudo->rhs) {ierr = VecDestroy(pseudo->rhs);CHKERRQ(ierr);} 208 if (ts->Ashell) {ierr = MatDestroy(ts->A); CHKERRQ(ierr);} 209 PetscFree(pseudo); 210 return 0; 211 } 212 213 214 /*------------------------------------------------------------*/ 215 /* 216 This matrix shell multiply where user provided Shell matrix 217 */ 218 219 int TSPseudoMatMult(Mat mat,Vec x,Vec y) 220 { 221 TS ts; 222 Scalar mdt,mone = -1.0; 223 int ierr; 224 225 MatShellGetContext(mat,(void **)&ts); 226 mdt = 1.0/ts->time_step; 227 228 /* apply user provided function */ 229 ierr = MatMult(ts->Ashell,x,y); CHKERRQ(ierr); 230 /* shift and scale by 1/dt - F */ 231 ierr = VecAXPBY(&mdt,&mone,x,y); CHKERRQ(ierr); 232 return 0; 233 } 234 235 /* 236 This defines the nonlinear equation that is to be solved with SNES 237 238 (U^{n+1} - U^{n})/dt - F(U^{n+1}) 239 */ 240 int TSPseudoFunction(SNES snes,Vec x,Vec y,void *ctx) 241 { 242 TS ts = (TS) ctx; 243 Scalar mdt = 1.0/ts->time_step,*unp1,*un,*Funp1; 244 int ierr,i,n; 245 246 /* apply user provided function */ 247 ierr = TSComputeRHSFunction(ts,ts->ptime,x,y); CHKERRQ(ierr); 248 /* compute (u^{n+1) - u^{n})/dt - F(u^{n+1}) */ 249 ierr = VecGetArray(ts->vec_sol,&un); CHKERRQ(ierr); 250 ierr = VecGetArray(x,&unp1); CHKERRQ(ierr); 251 ierr = VecGetArray(y,&Funp1); CHKERRQ(ierr); 252 ierr = VecGetLocalSize(x,&n); CHKERRQ(ierr); 253 for ( i=0; i<n; i++ ) { 254 Funp1[i] = mdt*(unp1[i] - un[i]) - Funp1[i]; 255 } 256 ierr = VecRestoreArray(ts->vec_sol,&un); 257 ierr = VecRestoreArray(x,&unp1); 258 ierr = VecRestoreArray(y,&Funp1); 259 260 return 0; 261 } 262 263 /* 264 This constructs the Jacobian needed for SNES 265 266 J = I/dt - J_{F} where J_{F} is the given Jacobian of F. 267 */ 268 int TSPseudoJacobian(SNES snes,Vec x,Mat *AA,Mat *BB,MatStructure *str,void *ctx) 269 { 270 TS ts = (TS) ctx; 271 int ierr; 272 Scalar mone = -1.0, mdt = 1.0/ts->time_step; 273 MatType mtype; 274 275 /* construct users Jacobian */ 276 if (ts->rhsjacobian) { 277 ierr = (*ts->rhsjacobian)(ts,ts->ptime,x,AA,BB,str,ts->jacP);CHKERRQ(ierr); 278 } 279 280 /* shift and scale Jacobian, if not a shell matrix */ 281 ierr = MatGetType(*AA,&mtype,PETSC_NULL); 282 if (mtype != MATSHELL) { 283 ierr = MatScale(&mone,*AA); CHKERRQ(ierr); 284 ierr = MatShift(&mdt,*AA); CHKERRQ(ierr); 285 } 286 ierr = MatGetType(*BB,&mtype,PETSC_NULL); 287 if (*BB != *AA && *str != SAME_PRECONDITIONER && mtype != MATSHELL) { 288 ierr = MatScale(&mone,*BB); CHKERRQ(ierr); 289 ierr = MatShift(&mdt,*BB); CHKERRQ(ierr); 290 } 291 292 return 0; 293 } 294 295 296 static int TSSetUp_Pseudo(TS ts) 297 { 298 TS_Pseudo *pseudo = (TS_Pseudo*) ts->data; 299 int ierr, M, m; 300 301 ierr = VecDuplicate(ts->vec_sol,&pseudo->update); CHKERRQ(ierr); 302 ierr = VecDuplicate(ts->vec_sol,&pseudo->func); CHKERRQ(ierr); 303 ierr = SNESSetFunction(ts->snes,pseudo->func,TSPseudoFunction,ts);CHKERRQ(ierr); 304 if (ts->Ashell) { /* construct new shell matrix */ 305 ierr = VecGetSize(ts->vec_sol,&M); CHKERRQ(ierr); 306 ierr = VecGetLocalSize(ts->vec_sol,&m); CHKERRQ(ierr); 307 ierr = MatCreateShell(ts->comm,m,M,M,M,ts,&ts->A); CHKERRQ(ierr); 308 ierr = MatShellSetOperation(ts->A,MATOP_MULT,(void*)TSPseudoMatMult);CHKERRQ(ierr); 309 } 310 ierr = SNESSetJacobian(ts->snes,ts->A,ts->B,TSPseudoJacobian,ts);CHKERRQ(ierr); 311 return 0; 312 } 313 /*------------------------------------------------------------*/ 314 315 int TSPseudoDefaultMonitor(TS ts, int step, double time,Vec v, void *ctx) 316 { 317 TS_Pseudo *pseudo = (TS_Pseudo*) ts->data; 318 319 PetscPrintf(ts->comm,"TS %d dt %g time %g fnorm %g\n",step,ts->time_step,time,pseudo->fnorm); 320 return 0; 321 } 322 323 static int TSSetFromOptions_Pseudo(TS ts) 324 { 325 int ierr,flg; 326 327 ierr = SNESSetFromOptions(ts->snes); CHKERRQ(ierr); 328 329 ierr = OptionsHasName(ts->prefix,"-ts_monitor",&flg); CHKERRQ(ierr); 330 if (flg) { 331 TSSetMonitor(ts,TSPseudoDefaultMonitor,0); 332 } 333 return 0; 334 } 335 336 static int TSPrintHelp_Pseudo(TS ts) 337 { 338 339 return 0; 340 } 341 342 static int TSView_Pseudo(PetscObject obj,Viewer viewer) 343 { 344 return 0; 345 } 346 347 /* ------------------------------------------------------------ */ 348 int TSCreate_Pseudo(TS ts ) 349 { 350 TS_Pseudo *pseudo; 351 int ierr; 352 MatType mtype; 353 354 ts->type = TS_PSEUDO; 355 ts->destroy = TSDestroy_Pseudo; 356 ts->printhelp = TSPrintHelp_Pseudo; 357 ts->view = TSView_Pseudo; 358 359 if (ts->problem_type == TS_LINEAR) { 360 SETERRQ(1,"TSCreate_Pseudo:Only for nonlinear problems"); 361 } 362 if (!ts->A) { 363 SETERRQ(1,"TSCreate_Pseudo:Must set Jacobian"); 364 } 365 ierr = MatGetType(ts->A,&mtype,PETSC_NULL); 366 if (mtype == MATSHELL) { 367 ts->Ashell = ts->A; 368 } 369 ts->setup = TSSetUp_Pseudo; 370 ts->step = TSStep_Pseudo; 371 ts->setfromoptions = TSSetFromOptions_Pseudo; 372 373 /* create the required nonlinear solver context */ 374 ierr = SNESCreate(ts->comm,SNES_NONLINEAR_EQUATIONS,&ts->snes);CHKERRQ(ierr); 375 376 pseudo = PetscNew(TS_Pseudo); CHKPTRQ(pseudo); 377 PetscMemzero(pseudo,sizeof(TS_Pseudo)); 378 ts->data = (void *) pseudo; 379 380 return 0; 381 } 382 383 384 385 386 387