Actual source code: beuler.c

  1: #define PETSCTS_DLL

  3: /*
  4:        Code for Timestepping with implicit backwards Euler.
  5: */
 6:  #include include/private/tsimpl.h

  8: typedef struct {
  9:   Vec  update;      /* work vector where new solution is formed */
 10:   Vec  func;        /* work vector where F(t[i],u[i]) is stored */
 11:   Vec  rhs;         /* work vector for RHS; vec_sol/dt */
 12: } TS_BEuler;

 14: /*------------------------------------------------------------------------------*/
 15: /* 
 16:    Set ts->A = ts->Arhs = 1/dt*Alhs - Arhs, used in KSPSolve() 
 17: */
 20: PetscErrorCode TSSetKSPOperators_BEuler(TS ts)
 21: {
 23:   PetscScalar    mdt = 1.0/ts->time_step;

 26:   if (!ts->A){
 27:     ts->A = ts->Arhs;
 28:     PetscObjectReference((PetscObject)ts->Arhs);
 29:   }

 31:   MatScale(ts->A,-1.0);
 32:   if (ts->Alhs){
 33:     MatAXPY(ts->A,mdt,ts->Alhs,DIFFERENT_NONZERO_PATTERN);
 34:   } else {
 35:     MatShift(ts->A,mdt);
 36:   }
 37:   return(0);
 38: }

 40: /*
 41:     Version for linear PDE where RHS does not depend on time. Has built a
 42:   single matrix that is to be used for all timesteps.
 43: */
 46: static PetscErrorCode TSStep_BEuler_Linear_Constant_Matrix(TS ts,PetscInt *steps,PetscReal *ptime)
 47: {
 48:   TS_BEuler      *beuler = (TS_BEuler*)ts->data;
 49:   Vec            sol = ts->vec_sol,update = beuler->update;
 50:   Vec            rhs = beuler->rhs;
 52:   PetscInt       i,max_steps = ts->max_steps,its;
 53:   PetscScalar    mdt = 1.0/ts->time_step;
 54:   KSP            ksp;

 57:   TSGetKSP(ts,&ksp);
 58:   *steps = -ts->steps;
 59:   TSMonitor(ts,ts->steps,ts->ptime,sol);

 61:   /* set initial guess to be previous solution */
 62:   VecCopy(sol,update);

 64:   for (i=0; i<max_steps; i++) {
 65:     /* set rhs = 1/dt*Alhs*sol */
 66:     if (ts->Alhs){
 67:       MatMult(ts->Alhs,sol,rhs);
 68:     } else {
 69:       VecCopy(sol,rhs);
 70:     }
 71:     VecScale(rhs,mdt);

 73:     ts->ptime += ts->time_step;
 74:     if (ts->ptime > ts->max_time) break;

 76:     /* solve (1/dt*Alhs - A)*update = rhs */
 77:     KSPSolve(ts->ksp,rhs,update);
 78:     KSPGetIterationNumber(ksp,&its);
 79:     ts->linear_its += its;
 80:     VecCopy(update,sol);
 81:     ts->steps++;
 82:     TSMonitor(ts,ts->steps,ts->ptime,sol);
 83:   }

 85:   *steps += ts->steps;
 86:   *ptime  = ts->ptime;
 87:   return(0);
 88: }

 90: /*
 91:       Version where matrix depends on time 
 92: */
 95: static PetscErrorCode TSStep_BEuler_Linear_Variable_Matrix(TS ts,PetscInt *steps,PetscReal *ptime)
 96: {
 97:   TS_BEuler      *beuler = (TS_BEuler*)ts->data;
 98:   Vec            sol = ts->vec_sol,update = beuler->update,rhs = beuler->rhs;
100:   PetscInt       i,max_steps = ts->max_steps,its;
101:   PetscScalar    mdt = 1.0/ts->time_step;
102:   PetscReal      t_mid;
103:   MatStructure   str;
104:   KSP            ksp;

107:   TSGetKSP(ts,&ksp);
108:   *steps = -ts->steps;
109:   TSMonitor(ts,ts->steps,ts->ptime,sol);

111:   /* set initial guess to be previous solution */
112:   VecCopy(sol,update);

114:   for (i=0; i<max_steps; i++) {
115:     /* set rhs = 1/dt*Alhs(t_mid)*sol */
116:     if (ts->Alhs){
117:       /* evaluate lhs matrix function at t_mid */
118:       t_mid = ts->ptime+ts->time_step/2.0;
119:       (*ts->ops->lhsmatrix)(ts,t_mid,&ts->Alhs,PETSC_NULL,&str,ts->jacPlhs);
120:       MatMult(ts->Alhs,sol,rhs);
121:     } else {
122:       VecCopy(sol,rhs);
123:     }
124:     VecScale(rhs,mdt);

126:     ts->ptime += ts->time_step;
127:     if (ts->ptime > ts->max_time) break;

129:     /* evaluate rhs matrix function at current ptime */
130:     (*ts->ops->rhsmatrix)(ts,ts->ptime,&ts->Arhs,&ts->B,&str,ts->jacP);

132:     /* set ts->A = ts->Arhs = 1/dt*Alhs - Arhs, used in KSPSolve() */
133:     TSSetKSPOperators_BEuler(ts);
134:     KSPSetOperators(ts->ksp,ts->A,ts->A,SAME_NONZERO_PATTERN);

136:     /* solve (1/dt*Alhs(t_mid) - A(t_n+1))*update = rhs */
137:     KSPSolve(ts->ksp,rhs,update);
138:     KSPGetIterationNumber(ksp,&its);
139:     ts->linear_its += its;
140:     VecCopy(update,sol);
141:     ts->steps++;
142:     TSMonitor(ts,ts->steps,ts->ptime,sol);
143:   }

145:   *steps += ts->steps;
146:   *ptime  = ts->ptime;
147:   return(0);
148: }
149: /*
150:     Version for nonlinear PDE.
151: */
154: static PetscErrorCode TSStep_BEuler_Nonlinear(TS ts,PetscInt *steps,PetscReal *ptime)
155: {
156:   Vec            sol = ts->vec_sol;
158:   PetscInt       i,max_steps = ts->max_steps,its,lits;
159:   TS_BEuler      *beuler = (TS_BEuler*)ts->data;
160: 
162:   *steps = -ts->steps;
163:   TSMonitor(ts,ts->steps,ts->ptime,sol);

165:   for (i=0; i<max_steps; i++) {
166:     ts->ptime += ts->time_step;
167:     if (ts->ptime > ts->max_time) break;
168:     VecCopy(sol,beuler->update);
169:     SNESSolve(ts->snes,PETSC_NULL,beuler->update);
170:     SNESGetNumberLinearIterations(ts->snes,&lits);
171:     SNESGetIterationNumber(ts->snes,&its);
172:     ts->nonlinear_its += its; ts->linear_its += lits;
173:     VecCopy(beuler->update,sol);
174:     ts->steps++;
175:     TSMonitor(ts,ts->steps,ts->ptime,sol);
176:   }

178:   *steps += ts->steps;
179:   *ptime  = ts->ptime;
180:   return(0);
181: }

183: /*------------------------------------------------------------*/
186: static PetscErrorCode TSDestroy_BEuler(TS ts)
187: {
188:   TS_BEuler      *beuler = (TS_BEuler*)ts->data;

192:   if (beuler->update) {VecDestroy(beuler->update);}
193:   if (beuler->func) {VecDestroy(beuler->func);}
194:   if (beuler->rhs) {VecDestroy(beuler->rhs);}
195:   PetscFree(beuler);
196:   return(0);
197: }

199: /* 
200:     This defines the nonlinear equation that is to be solved with SNES
201:       1/dt* (U^{n+1} - U^{n}) - F(U^{n+1}) 
202: */
205: PetscErrorCode TSBEulerFunction(SNES snes,Vec x,Vec y,void *ctx)
206: {
207:   TS             ts = (TS) ctx;
208:   PetscScalar    mdt = 1.0/ts->time_step,*unp1,*un,*Funp1;
210:   PetscInt       i,n;

213:   /* apply user-provided function */
214:   TSComputeRHSFunction(ts,ts->ptime,x,y);
215:   VecGetArray(ts->vec_sol,&un);
216:   VecGetArray(x,&unp1);
217:   VecGetArray(y,&Funp1);
218:   VecGetLocalSize(x,&n);
219:   for (i=0; i<n; i++) {
220:     Funp1[i] = mdt*(unp1[i] - un[i]) - Funp1[i];
221:   }
222:   VecRestoreArray(ts->vec_sol,&un);
223:   VecRestoreArray(x,&unp1);
224:   VecRestoreArray(y,&Funp1);
225:   return(0);
226: }

228: /*
229:    This constructs the Jacobian needed for SNES 
230:      J = I/dt - J_{F}   where J_{F} is the given Jacobian of F at t_{n+1}.
231:      x  - input vector
232:      AA - Jacobian matrix 
233:      BB - preconditioner matrix, usually the same as AA
234: */
237: PetscErrorCode TSBEulerJacobian(SNES snes,Vec x,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
238: {
239:   TS             ts = (TS) ctx;

243:   /* construct user's Jacobian */
244:   TSComputeRHSJacobian(ts,ts->ptime,x,AA,BB,str);

246:   /* shift and scale Jacobian */
247:   /* this test is a undesirable hack, we assume that if it is MATMFFD then it is
248:      obtained from -snes_mf_operator and there is computed directly from the 
249:      FormFunction() SNES is given and therefor does not need to be shifted/scaled
250:      BUT maybe it could be MATMFFD and does require shift in some other case??? */
251:   TSSetKSPOperators_BEuler(ts);
252:   return(0);
253: }

255: /* ------------------------------------------------------------*/
258: static PetscErrorCode TSSetUp_BEuler_Linear_Constant_Matrix(TS ts)
259: {
260:   TS_BEuler      *beuler = (TS_BEuler*)ts->data;

264:   KSPSetFromOptions(ts->ksp);
265:   VecDuplicate(ts->vec_sol,&beuler->update);
266:   VecDuplicate(ts->vec_sol,&beuler->rhs);
267: 
268:   /* build linear system to be solved - should move into TSStep() if dt changes! */
269:   /* Set ts->A = ts->Arhs = 1/dt*Alhs - Arhs, used in KSPSolve() */
270:   TSSetKSPOperators_BEuler(ts);
271:   KSPSetOperators(ts->ksp,ts->A,ts->A,SAME_NONZERO_PATTERN);
272:   return(0);
273: }

277: static PetscErrorCode TSSetUp_BEuler_Linear_Variable_Matrix(TS ts)
278: {
279:   TS_BEuler      *beuler = (TS_BEuler*)ts->data;

283:   KSPSetFromOptions(ts->ksp);
284:   VecDuplicate(ts->vec_sol,&beuler->update);
285:   VecDuplicate(ts->vec_sol,&beuler->rhs);
286:   return(0);
287: }

291: static PetscErrorCode TSSetUp_BEuler_Nonlinear(TS ts)
292: {
293:   TS_BEuler      *beuler = (TS_BEuler*)ts->data;

297:   VecDuplicate(ts->vec_sol,&beuler->update);
298:   VecDuplicate(ts->vec_sol,&beuler->func);
299:   SNESSetFunction(ts->snes,beuler->func,TSBEulerFunction,ts);
300:   SNESSetJacobian(ts->snes,ts->Arhs,ts->B,TSBEulerJacobian,ts);
301:   return(0);
302: }
303: /*------------------------------------------------------------*/

307: static PetscErrorCode TSSetFromOptions_BEuler_Linear(TS ts)
308: {
310:   return(0);
311: }

315: static PetscErrorCode TSSetFromOptions_BEuler_Nonlinear(TS ts)
316: {
318:   return(0);
319: }

323: static PetscErrorCode TSView_BEuler(TS ts,PetscViewer viewer)
324: {
326:   return(0);
327: }

329: /* ------------------------------------------------------------ */
330: /*MC
331:       TS_BEULER - ODE solver using the implicit backward Euler method

333:   Level: beginner

335: .seealso:  TSCreate(), TS, TSSetType(), TS_EULER

337: M*/
341: PetscErrorCode  TSCreate_BEuler(TS ts)
342: {
343:   TS_BEuler      *beuler;

347:   ts->ops->destroy = TSDestroy_BEuler;
348:   ts->ops->view    = TSView_BEuler;

350:   if (ts->problem_type == TS_LINEAR) {
351:     if (!ts->Arhs) {
352:       SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set rhs matrix for linear problem");
353:     }
354:     if (!ts->ops->rhsmatrix) {
355:       ts->ops->setup  = TSSetUp_BEuler_Linear_Constant_Matrix;
356:       ts->ops->step   = TSStep_BEuler_Linear_Constant_Matrix;
357:     } else {
358:       ts->ops->setup  = TSSetUp_BEuler_Linear_Variable_Matrix;
359:       ts->ops->step   = TSStep_BEuler_Linear_Variable_Matrix;
360:     }
361:     ts->ops->setfromoptions  = TSSetFromOptions_BEuler_Linear;
362:     KSPCreate(ts->comm,&ts->ksp);
363:     KSPSetInitialGuessNonzero(ts->ksp,PETSC_TRUE);
364:   } else if (ts->problem_type == TS_NONLINEAR) {
365:     ts->ops->setup           = TSSetUp_BEuler_Nonlinear;
366:     ts->ops->step            = TSStep_BEuler_Nonlinear;
367:     ts->ops->setfromoptions  = TSSetFromOptions_BEuler_Nonlinear;
368:     SNESCreate(ts->comm,&ts->snes);
369:   } else SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"No such problem");

371:   PetscNew(TS_BEuler,&beuler);
372:   PetscLogObjectMemory(ts,sizeof(TS_BEuler));
373:   ts->data = (void*)beuler;

375:   return(0);
376: }