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*); 20 void *dtctx; 21 double initial_fnorm,fnorm; /* original and current norm of F(u) */ 22 23 } TS_PosIndep; 24 25 /* ------------------------------------------------------------------------------*/ 26 /*@C 27 TSPseudoDefaultPosIndTimeStep - Default code to compute 28 pseudo timestepping. Use with TSPseudoSetPosIndTimeStep(). 29 30 Input Parameters: 31 . ts - the time step context 32 . dtctx - 33 34 Output Parameter: 35 . newdt - the time step to use for the next step 36 37 @*/ 38 int TSPseudoDefaultPosIndTimeStep(TS ts,double* newdt,void* dtctx) 39 { 40 TS_PosIndep *posindep = (TS_PosIndep*) ts->data; 41 int ierr; 42 43 ierr = TSComputeRHSFunction(ts,ts->ptime,ts->vec_sol,posindep->func);CHKERRQ(ierr); 44 ierr = VecNorm(posindep->func,NORM_2,&posindep->fnorm); CHKERRQ(ierr); 45 if (posindep->initial_fnorm == 0.0) { 46 /* first time through so compute initial function norm */ 47 posindep->initial_fnorm = posindep->fnorm; 48 } 49 if (posindep->fnorm == 0.0) *newdt = 1.e12*1.1*ts->time_step; 50 else *newdt = 1.1*ts->time_step*posindep->initial_fnorm/posindep->fnorm; 51 return 0; 52 } 53 54 /*@ 55 TSPseudoSetPosIndTimeStep - Sets the user routine to be 56 called at each pseudo-time-step to update the time-step. 57 58 Input Parameters: 59 . ts - time step context 60 . dt - function to compute timestep 61 . ctx - [optional] context required by function 62 63 @*/ 64 int TSPseudoSetPosIndTimeStep(TS ts,int (*dt)(TS,double*,void*),void* ctx) 65 { 66 TS_PosIndep *posindep; 67 68 PetscValidHeaderSpecific(ts,TS_COOKIE); 69 if (ts->type != TS_PSEUDO_POSIND) return 0; 70 posindep = (TS_PosIndep*) ts->data; 71 posindep->dt = dt; 72 posindep->dtctx = ctx; 73 return 0; 74 } 75 76 /* 77 78 */ 79 static int TSStep_PosIndep(TS ts,int *steps,double *time) 80 { 81 Vec sol = ts->vec_sol; 82 int ierr,i,max_steps = ts->max_steps,its; 83 TS_PosIndep *posindep = (TS_PosIndep*) ts->data; 84 85 *steps = -ts->steps; 86 87 ierr = VecCopy(sol,posindep->update); CHKERRQ(ierr); 88 for ( i=0; i<max_steps && ts->ptime < ts->max_time; i++ ) { 89 ierr = (*posindep->dt)(ts,&ts->time_step, posindep->dtctx); CHKERRQ(ierr); 90 ts->ptime += ts->time_step; 91 ierr = SNESSolve(ts->snes,posindep->update,&its); CHKERRQ(ierr); 92 ierr = VecCopy(posindep->update,sol); CHKERRQ(ierr); 93 ts->steps++; 94 ierr = TSMonitor(ts,ts->steps,ts->ptime,sol);CHKERRQ(ierr); 95 } 96 97 *steps += ts->steps; 98 *time = ts->ptime; 99 return 0; 100 } 101 102 /*------------------------------------------------------------*/ 103 static int TSDestroy_PosIndep(PetscObject obj ) 104 { 105 TS ts = (TS) obj; 106 TS_PosIndep *posindep = (TS_PosIndep*) ts->data; 107 int ierr; 108 109 ierr = VecDestroy(posindep->update); CHKERRQ(ierr); 110 if (posindep->func) {ierr = VecDestroy(posindep->func);CHKERRQ(ierr);} 111 if (posindep->rhs) {ierr = VecDestroy(posindep->rhs);CHKERRQ(ierr);} 112 if (ts->Ashell) {ierr = MatDestroy(ts->A); CHKERRQ(ierr);} 113 PetscFree(posindep); 114 return 0; 115 } 116 117 118 /*------------------------------------------------------------*/ 119 /* 120 This matrix shell multiply where user provided Shell matrix 121 */ 122 123 int TSPosIndepMatMult(Mat mat,Vec x,Vec y) 124 { 125 TS ts; 126 Scalar mdt,mone = -1.0; 127 int ierr; 128 129 MatShellGetContext(mat,(void **)&ts); 130 mdt = 1.0/ts->time_step; 131 132 /* apply user provided function */ 133 ierr = MatMult(ts->Ashell,x,y); CHKERRQ(ierr); 134 /* shift and scale by 1/dt - F */ 135 ierr = VecAXPBY(&mdt,&mone,x,y); CHKERRQ(ierr); 136 return 0; 137 } 138 139 /* 140 This defines the nonlinear equation that is to be solved with SNES 141 142 (U^{n+1} - U^{n})/dt - F(U^{n+1}) 143 */ 144 int TSPosIndepFunction(SNES snes,Vec x,Vec y,void *ctx) 145 { 146 TS ts = (TS) ctx; 147 Scalar mdt = 1.0/ts->time_step,*unp1,*un,*Funp1; 148 int ierr,i,n; 149 150 /* apply user provided function */ 151 ierr = TSComputeRHSFunction(ts,ts->ptime,x,y); CHKERRQ(ierr); 152 /* (u^{n+1) - u^{n})/dt - F(u^{n+1}) */ 153 ierr = VecGetArray(ts->vec_sol,&un); CHKERRQ(ierr); 154 ierr = VecGetArray(x,&unp1); CHKERRQ(ierr); 155 ierr = VecGetArray(y,&Funp1); CHKERRQ(ierr); 156 ierr = VecGetLocalSize(x,&n); CHKERRQ(ierr); 157 for ( i=0; i<n; i++ ) { 158 Funp1[i] = mdt*(unp1[i] - un[i]) - Funp1[i]; 159 } 160 ierr = VecRestoreArray(ts->vec_sol,&un); 161 ierr = VecRestoreArray(x,&unp1); 162 ierr = VecRestoreArray(y,&Funp1); 163 164 return 0; 165 } 166 167 /* 168 This constructs the Jacobian needed for SNES 169 170 J = I/dt - J_{F} where J_{F} is the given Jacobian of F. 171 */ 172 int TSPosIndepJacobian(SNES snes,Vec x,Mat *AA,Mat *BB,MatStructure *str,void *ctx) 173 { 174 TS ts = (TS) ctx; 175 int ierr; 176 Scalar mone = -1.0, mdt = 1.0/ts->time_step; 177 MatType mtype; 178 179 /* construct users Jacobian */ 180 if (ts->rhsjacobian) { 181 ierr = (*ts->rhsjacobian)(ts,ts->ptime,x,AA,BB,str,ts->jacP);CHKERRQ(ierr); 182 } 183 184 /* shift and scale Jacobian, if not a shell matrix */ 185 ierr = MatGetType(*AA,&mtype,PETSC_NULL); 186 if (mtype != MATSHELL) { 187 ierr = MatScale(&mone,*AA); CHKERRQ(ierr); 188 ierr = MatShift(&mdt,*AA); CHKERRQ(ierr); 189 } 190 ierr = MatGetType(*BB,&mtype,PETSC_NULL); 191 if (*BB != *AA && *str != SAME_PRECONDITIONER && mtype != MATSHELL) { 192 ierr = MatScale(&mone,*BB); CHKERRQ(ierr); 193 ierr = MatShift(&mdt,*BB); CHKERRQ(ierr); 194 } 195 196 return 0; 197 } 198 199 200 static int TSSetUp_PosIndep(TS ts) 201 { 202 TS_PosIndep *posindep = (TS_PosIndep*) ts->data; 203 int ierr, M, m; 204 205 ierr = VecDuplicate(ts->vec_sol,&posindep->update); CHKERRQ(ierr); 206 ierr = VecDuplicate(ts->vec_sol,&posindep->func); CHKERRQ(ierr); 207 ierr = SNESSetFunction(ts->snes,posindep->func,TSPosIndepFunction,ts);CHKERRQ(ierr); 208 if (ts->Ashell) { /* construct new shell matrix */ 209 ierr = VecGetSize(ts->vec_sol,&M); CHKERRQ(ierr); 210 ierr = VecGetLocalSize(ts->vec_sol,&m); CHKERRQ(ierr); 211 ierr = MatCreateShell(ts->comm,m,M,M,M,ts,&ts->A); CHKERRQ(ierr); 212 ierr = MatShellSetOperation(ts->A,MAT_MULT,(void*)TSPosIndepMatMult);CHKERRQ(ierr); 213 } 214 ierr = SNESSetJacobian(ts->snes,ts->A,ts->B,TSPosIndepJacobian,ts);CHKERRQ(ierr); 215 return 0; 216 } 217 /*------------------------------------------------------------*/ 218 219 int TSPseudoDefaultMonitor(TS ts, int step, double time,Vec v, void *ctx) 220 { 221 TS_PosIndep *posindep = (TS_PosIndep*) ts->data; 222 223 PetscPrintf(ts->comm,"TS %d dt %g time %g fnorm %g\n",step,ts->time_step,time, 224 posindep->fnorm); 225 return 0; 226 } 227 228 229 static int TSSetFromOptions_PosIndep(TS ts) 230 { 231 int ierr,flg; 232 233 ierr = SNESSetFromOptions(ts->snes); CHKERRQ(ierr); 234 235 ierr = OptionsHasName(ts->prefix,"-ts_monitor",&flg); CHKERRQ(ierr); 236 if (flg) { 237 TSSetMonitor(ts,TSPseudoDefaultMonitor,0); 238 } 239 return 0; 240 } 241 242 static int TSPrintHelp_PosIndep(TS ts) 243 { 244 245 return 0; 246 } 247 248 static int TSView_PosIndep(PetscObject obj,Viewer viewer) 249 { 250 return 0; 251 } 252 253 /* ------------------------------------------------------------ */ 254 int TSCreate_PosIndep(TS ts ) 255 { 256 TS_PosIndep *posindep; 257 int ierr; 258 MatType mtype; 259 260 ts->type = TS_PSEUDO_POSIND; 261 ts->destroy = TSDestroy_PosIndep; 262 ts->printhelp = TSPrintHelp_PosIndep; 263 ts->view = TSView_PosIndep; 264 265 if (ts->problem_type == TS_LINEAR) { 266 SETERRQ(1,"TSCreate_PosIndep:Only for nonlinear problems"); 267 } 268 if (!ts->A) { 269 SETERRQ(1,"TSCreate_PosIndep:Must set Jacobian"); 270 } 271 ierr = MatGetType(ts->A,&mtype,PETSC_NULL); 272 if (mtype == MATSHELL) { 273 ts->Ashell = ts->A; 274 } 275 ts->setup = TSSetUp_PosIndep; 276 ts->step = TSStep_PosIndep; 277 ts->setfromoptions = TSSetFromOptions_PosIndep; 278 ierr = SNESCreate(ts->comm,SNES_NONLINEAR_EQUATIONS,&ts->snes);CHKERRQ(ierr); 279 280 posindep = PetscNew(TS_PosIndep); CHKPTRQ(posindep); 281 PetscMemzero(posindep,sizeof(TS_PosIndep)); 282 ts->data = (void *) posindep; 283 284 return 0; 285 } 286 287 288 289 290 291