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