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: }