1 static char help[] = "Chemo-taxis Problems from Mathematical Biology.\n";
2
3 /*
4 Page 18, Chemo-taxis Problems from Mathematical Biology
5
6 rho_t =
7 c_t =
8
9 Further discussion on Page 134 and in 2d on Page 409
10 */
11
12 /*
13
14 Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
15 Include "petscts.h" so that we can use SNES solvers. Note that this
16 file automatically includes:
17 petscsys.h - base PETSc routines petscvec.h - vectors
18 petscmat.h - matrices
19 petscis.h - index sets petscksp.h - Krylov subspace methods
20 petscviewer.h - viewers petscpc.h - preconditioners
21 petscksp.h - linear solvers
22 */
23
24 #include <petscdm.h>
25 #include <petscdmda.h>
26 #include <petscts.h>
27
28 typedef struct {
29 PetscScalar rho, c;
30 } Field;
31
32 typedef struct {
33 PetscScalar epsilon, delta, alpha, beta, gamma, kappa, lambda, mu, cstar;
34 PetscBool upwind;
35 } AppCtx;
36
37 /*
38 User-defined routines
39 */
40 extern PetscErrorCode IFunction(TS, PetscReal, Vec, Vec, Vec, void *), InitialConditions(DM, Vec);
41
main(int argc,char ** argv)42 int main(int argc, char **argv)
43 {
44 TS ts; /* nonlinear solver */
45 Vec U; /* solution, residual vectors */
46 DM da;
47 AppCtx appctx;
48
49 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50 Initialize program
51 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
52 PetscFunctionBeginUser;
53 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
54
55 appctx.epsilon = 1.0e-3;
56 appctx.delta = 1.0;
57 appctx.alpha = 10.0;
58 appctx.beta = 4.0;
59 appctx.gamma = 1.0;
60 appctx.kappa = .75;
61 appctx.lambda = 1.0;
62 appctx.mu = 100.;
63 appctx.cstar = .2;
64 appctx.upwind = PETSC_TRUE;
65
66 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-delta", &appctx.delta, NULL));
67 PetscCall(PetscOptionsGetBool(NULL, NULL, "-upwind", &appctx.upwind, NULL));
68
69 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70 Create distributed array (DMDA) to manage parallel grid and vectors
71 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
72 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 8, 2, 1, NULL, &da));
73 PetscCall(DMSetFromOptions(da));
74 PetscCall(DMSetUp(da));
75 PetscCall(DMDASetFieldName(da, 0, "rho"));
76 PetscCall(DMDASetFieldName(da, 1, "c"));
77
78 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79 Extract global vectors from DMDA; then duplicate for remaining
80 vectors that are the same types
81 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
82 PetscCall(DMCreateGlobalVector(da, &U));
83
84 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85 Create timestepping solver context
86 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
88 PetscCall(TSSetType(ts, TSROSW));
89 PetscCall(TSSetDM(ts, da));
90 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
91 PetscCall(TSSetIFunction(ts, NULL, IFunction, &appctx));
92
93 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94 Set initial conditions
95 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96 PetscCall(InitialConditions(da, U));
97 PetscCall(TSSetSolution(ts, U));
98
99 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100 Set solver options
101 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102 PetscCall(TSSetTimeStep(ts, .0001));
103 PetscCall(TSSetMaxTime(ts, 1.0));
104 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
105 PetscCall(TSSetFromOptions(ts));
106
107 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108 Solve nonlinear system
109 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110 PetscCall(TSSolve(ts, U));
111
112 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113 Free work space. All PETSc objects should be destroyed when they
114 are no longer needed.
115 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116 PetscCall(VecDestroy(&U));
117 PetscCall(TSDestroy(&ts));
118 PetscCall(DMDestroy(&da));
119
120 PetscCall(PetscFinalize());
121 return 0;
122 }
123 /* ------------------------------------------------------------------- */
124 /*
125 IFunction - Evaluates nonlinear function, F(U).
126
127 Input Parameters:
128 . ts - the TS context
129 . U - input vector
130 . ptr - optional user-defined context, as set by SNESSetFunction()
131
132 Output Parameter:
133 . F - function vector
134 */
IFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void * ptr)135 PetscErrorCode IFunction(TS ts, PetscReal ftime, Vec U, Vec Udot, Vec F, void *ptr)
136 {
137 AppCtx *appctx = (AppCtx *)ptr;
138 DM da;
139 PetscInt i, Mx, xs, xm;
140 PetscReal hx, sx;
141 PetscScalar rho, c, rhoxx, cxx, cx, rhox, kcxrhox;
142 Field *u, *f, *udot;
143 Vec localU;
144
145 PetscFunctionBegin;
146 PetscCall(TSGetDM(ts, &da));
147 PetscCall(DMGetLocalVector(da, &localU));
148 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
149
150 hx = 1.0 / (PetscReal)(Mx - 1);
151 sx = 1.0 / (hx * hx);
152
153 /*
154 Scatter ghost points to local vector,using the 2-step process
155 DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
156 By placing code between these two statements, computations can be
157 done while messages are in transition.
158 */
159 PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
160 PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));
161
162 /*
163 Get pointers to vector data
164 */
165 PetscCall(DMDAVecGetArrayRead(da, localU, &u));
166 PetscCall(DMDAVecGetArrayRead(da, Udot, &udot));
167 PetscCall(DMDAVecGetArrayWrite(da, F, &f));
168
169 /*
170 Get local grid boundaries
171 */
172 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
173
174 if (!xs) {
175 f[0].rho = udot[0].rho; /* u[0].rho - 0.0; */
176 f[0].c = udot[0].c; /* u[0].c - 1.0; */
177 xs++;
178 xm--;
179 }
180 if (xs + xm == Mx) {
181 f[Mx - 1].rho = udot[Mx - 1].rho; /* u[Mx-1].rho - 1.0; */
182 f[Mx - 1].c = udot[Mx - 1].c; /* u[Mx-1].c - 0.0; */
183 xm--;
184 }
185
186 /*
187 Compute function over the locally owned part of the grid
188 */
189 for (i = xs; i < xs + xm; i++) {
190 rho = u[i].rho;
191 rhoxx = (-2.0 * rho + u[i - 1].rho + u[i + 1].rho) * sx;
192 c = u[i].c;
193 cxx = (-2.0 * c + u[i - 1].c + u[i + 1].c) * sx;
194
195 if (!appctx->upwind) {
196 rhox = .5 * (u[i + 1].rho - u[i - 1].rho) / hx;
197 cx = .5 * (u[i + 1].c - u[i - 1].c) / hx;
198 kcxrhox = appctx->kappa * (cxx * rho + cx * rhox);
199 } else {
200 kcxrhox = appctx->kappa * ((u[i + 1].c - u[i].c) * u[i + 1].rho - (u[i].c - u[i - 1].c) * u[i].rho) * sx;
201 }
202
203 f[i].rho = udot[i].rho - appctx->epsilon * rhoxx + kcxrhox - appctx->mu * PetscAbsScalar(rho) * (1.0 - rho) * PetscMax(0, PetscRealPart(c - appctx->cstar)) + appctx->beta * rho;
204 f[i].c = udot[i].c - appctx->delta * cxx + appctx->lambda * c + appctx->alpha * rho * c / (appctx->gamma + c);
205 }
206
207 /*
208 Restore vectors
209 */
210 PetscCall(DMDAVecRestoreArrayRead(da, localU, &u));
211 PetscCall(DMDAVecRestoreArrayRead(da, Udot, &udot));
212 PetscCall(DMDAVecRestoreArrayWrite(da, F, &f));
213 PetscCall(DMRestoreLocalVector(da, &localU));
214 PetscFunctionReturn(PETSC_SUCCESS);
215 }
216
217 /* ------------------------------------------------------------------- */
InitialConditions(DM da,Vec U)218 PetscErrorCode InitialConditions(DM da, Vec U)
219 {
220 PetscInt i, xs, xm, Mx;
221 Field *u;
222 PetscReal hx, x;
223
224 PetscFunctionBegin;
225 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
226
227 hx = 1.0 / (PetscReal)(Mx - 1);
228
229 /*
230 Get pointers to vector data
231 */
232 PetscCall(DMDAVecGetArrayWrite(da, U, &u));
233
234 /*
235 Get local grid boundaries
236 */
237 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
238
239 /*
240 Compute function over the locally owned part of the grid
241 */
242 for (i = xs; i < xs + xm; i++) {
243 x = i * hx;
244 if (i < Mx - 1) u[i].rho = 0.0;
245 else u[i].rho = 1.0;
246 u[i].c = PetscCosReal(.5 * PETSC_PI * x);
247 }
248
249 /*
250 Restore vectors
251 */
252 PetscCall(DMDAVecRestoreArrayWrite(da, U, &u));
253 PetscFunctionReturn(PETSC_SUCCESS);
254 }
255
256 /*TEST
257
258 test:
259 args: -pc_type lu -da_refine 2 -ts_view -ts_monitor -ts_max_time 1
260 requires: double
261
262 test:
263 suffix: 2
264 args: -pc_type lu -da_refine 2 -ts_view -ts_monitor_draw_solution -ts_monitor -ts_max_time 1
265 requires: x double
266 output_file: output/ex4_1.out
267
268 TEST*/
269