1c4762a1bSJed Brown static char help[] = "Nonlinear Radiative Transport PDE with multigrid in 2d.\n\
2c4762a1bSJed Brown Uses 2-dimensional distributed arrays.\n\
3c4762a1bSJed Brown A 2-dim simplified Radiative Transport test problem is used, with analytic Jacobian. \n\
4c4762a1bSJed Brown \n\
5c4762a1bSJed Brown Solves the linear systems via multilevel methods \n\
6c4762a1bSJed Brown \n\
7c4762a1bSJed Brown The command line\n\
8c4762a1bSJed Brown options are:\n\
9c4762a1bSJed Brown -tleft <tl>, where <tl> indicates the left Diriclet BC \n\
10c4762a1bSJed Brown -tright <tr>, where <tr> indicates the right Diriclet BC \n\
11c4762a1bSJed Brown -beta <beta>, where <beta> indicates the exponent in T \n\n";
12c4762a1bSJed Brown
13c4762a1bSJed Brown /*
14c4762a1bSJed Brown
15c4762a1bSJed Brown This example models the partial differential equation
16c4762a1bSJed Brown
17c4762a1bSJed Brown - Div(alpha* T^beta (GRAD T)) = 0.
18c4762a1bSJed Brown
19c4762a1bSJed Brown where beta = 2.5 and alpha = 1.0
20c4762a1bSJed Brown
21c4762a1bSJed Brown BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0.
22c4762a1bSJed Brown
23c4762a1bSJed Brown in the unit square, which is uniformly discretized in each of x and
24c4762a1bSJed Brown y in this simple encoding. The degrees of freedom are cell centered.
25c4762a1bSJed Brown
26c4762a1bSJed Brown A finite volume approximation with the usual 5-point stencil
27c4762a1bSJed Brown is used to discretize the boundary value problem to obtain a
28c4762a1bSJed Brown nonlinear system of equations.
29c4762a1bSJed Brown
30c4762a1bSJed Brown This code was contributed by David Keyes
31c4762a1bSJed Brown
32c4762a1bSJed Brown */
33c4762a1bSJed Brown
34c4762a1bSJed Brown #include <petscsnes.h>
35c4762a1bSJed Brown #include <petscdm.h>
36c4762a1bSJed Brown #include <petscdmda.h>
37c4762a1bSJed Brown
38c4762a1bSJed Brown /* User-defined application context */
39c4762a1bSJed Brown
40c4762a1bSJed Brown typedef struct {
41c4762a1bSJed Brown PetscReal tleft, tright; /* Dirichlet boundary conditions */
42c4762a1bSJed Brown PetscReal beta, bm1, coef; /* nonlinear diffusivity parameterizations */
43c4762a1bSJed Brown } AppCtx;
44c4762a1bSJed Brown
45c4762a1bSJed Brown #define POWFLOP 5 /* assume a pow() takes five flops */
46c4762a1bSJed Brown
47c4762a1bSJed Brown extern PetscErrorCode FormInitialGuess(SNES, Vec, void *);
48c4762a1bSJed Brown extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
49c4762a1bSJed Brown extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
50c4762a1bSJed Brown
main(int argc,char ** argv)51d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv)
52d71ae5a4SJacob Faibussowitsch {
53c4762a1bSJed Brown SNES snes;
54c4762a1bSJed Brown AppCtx user;
55c4762a1bSJed Brown PetscInt its, lits;
56c4762a1bSJed Brown PetscReal litspit;
57c4762a1bSJed Brown DM da;
58c4762a1bSJed Brown
59327415f7SBarry Smith PetscFunctionBeginUser;
609566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &argv, NULL, help));
61c4762a1bSJed Brown
62c4762a1bSJed Brown /* set problem parameters */
63c4762a1bSJed Brown user.tleft = 1.0;
64c4762a1bSJed Brown user.tright = 0.1;
65c4762a1bSJed Brown user.beta = 2.5;
669566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetReal(NULL, NULL, "-tleft", &user.tleft, NULL));
679566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetReal(NULL, NULL, "-tright", &user.tright, NULL));
689566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetReal(NULL, NULL, "-beta", &user.beta, NULL));
69c4762a1bSJed Brown user.bm1 = user.beta - 1.0;
70c4762a1bSJed Brown user.coef = user.beta / 2.0;
71c4762a1bSJed Brown
72c4762a1bSJed Brown /*
73c4762a1bSJed Brown Create the multilevel DM data structure
74c4762a1bSJed Brown */
759566063dSJacob Faibussowitsch PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
76c4762a1bSJed Brown
77c4762a1bSJed Brown /*
78c4762a1bSJed Brown Set the DMDA (grid structure) for the grids.
79c4762a1bSJed Brown */
809566063dSJacob Faibussowitsch PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 5, 5, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, &da));
819566063dSJacob Faibussowitsch PetscCall(DMSetFromOptions(da));
829566063dSJacob Faibussowitsch PetscCall(DMSetUp(da));
839566063dSJacob Faibussowitsch PetscCall(DMSetApplicationContext(da, &user));
849566063dSJacob Faibussowitsch PetscCall(SNESSetDM(snes, (DM)da));
85c4762a1bSJed Brown
86c4762a1bSJed Brown /*
87c4762a1bSJed Brown Create the nonlinear solver, and tell it the functions to use
88c4762a1bSJed Brown */
899566063dSJacob Faibussowitsch PetscCall(SNESSetFunction(snes, NULL, FormFunction, &user));
909566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, NULL, NULL, FormJacobian, &user));
919566063dSJacob Faibussowitsch PetscCall(SNESSetFromOptions(snes));
929566063dSJacob Faibussowitsch PetscCall(SNESSetComputeInitialGuess(snes, FormInitialGuess, NULL));
93c4762a1bSJed Brown
949566063dSJacob Faibussowitsch PetscCall(SNESSolve(snes, NULL, NULL));
959566063dSJacob Faibussowitsch PetscCall(SNESGetIterationNumber(snes, &its));
969566063dSJacob Faibussowitsch PetscCall(SNESGetLinearSolveIterations(snes, &lits));
97c4762a1bSJed Brown litspit = ((PetscReal)lits) / ((PetscReal)its);
9863a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
9963a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of Linear iterations = %" PetscInt_FMT "\n", lits));
1009566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Average Linear its / SNES = %e\n", (double)litspit));
101c4762a1bSJed Brown
1029566063dSJacob Faibussowitsch PetscCall(DMDestroy(&da));
1039566063dSJacob Faibussowitsch PetscCall(SNESDestroy(&snes));
1049566063dSJacob Faibussowitsch PetscCall(PetscFinalize());
105b122ec5aSJacob Faibussowitsch return 0;
106c4762a1bSJed Brown }
107c4762a1bSJed Brown /* -------------------- Form initial approximation ----------------- */
FormInitialGuess(SNES snes,Vec X,PetscCtx ctx)108*2a8381b2SBarry Smith PetscErrorCode FormInitialGuess(SNES snes, Vec X, PetscCtx ctx)
109d71ae5a4SJacob Faibussowitsch {
110c4762a1bSJed Brown AppCtx *user;
111c4762a1bSJed Brown PetscInt i, j, xs, ys, xm, ym;
112c4762a1bSJed Brown PetscReal tleft;
113c4762a1bSJed Brown PetscScalar **x;
114c4762a1bSJed Brown DM da;
115c4762a1bSJed Brown
116c4762a1bSJed Brown PetscFunctionBeginUser;
1179566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &da));
1189566063dSJacob Faibussowitsch PetscCall(DMGetApplicationContext(da, &user));
119c4762a1bSJed Brown tleft = user->tleft;
120c4762a1bSJed Brown /* Get ghost points */
1219566063dSJacob Faibussowitsch PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0));
1229566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArray(da, X, &x));
123c4762a1bSJed Brown
124c4762a1bSJed Brown /* Compute initial guess */
125c4762a1bSJed Brown for (j = ys; j < ys + ym; j++) {
126ad540459SPierre Jolivet for (i = xs; i < xs + xm; i++) x[j][i] = tleft;
127c4762a1bSJed Brown }
1289566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArray(da, X, &x));
1293ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
130c4762a1bSJed Brown }
131c4762a1bSJed Brown /* -------------------- Evaluate Function F(x) --------------------- */
FormFunction(SNES snes,Vec X,Vec F,void * ptr)132d71ae5a4SJacob Faibussowitsch PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr)
133d71ae5a4SJacob Faibussowitsch {
134c4762a1bSJed Brown AppCtx *user = (AppCtx *)ptr;
135c4762a1bSJed Brown PetscInt i, j, mx, my, xs, ys, xm, ym;
136c4762a1bSJed Brown PetscScalar zero = 0.0, one = 1.0;
137c4762a1bSJed Brown PetscScalar hx, hy, hxdhy, hydhx;
138c4762a1bSJed Brown PetscScalar t0, tn, ts, te, tw, an, as, ae, aw, dn, ds, de, dw, fn = 0.0, fs = 0.0, fe = 0.0, fw = 0.0;
139c4762a1bSJed Brown PetscScalar tleft, tright, beta;
140c4762a1bSJed Brown PetscScalar **x, **f;
141c4762a1bSJed Brown Vec localX;
142c4762a1bSJed Brown DM da;
143c4762a1bSJed Brown
144c4762a1bSJed Brown PetscFunctionBeginUser;
1459566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &da));
1469566063dSJacob Faibussowitsch PetscCall(DMGetLocalVector(da, &localX));
1479566063dSJacob Faibussowitsch PetscCall(DMDAGetInfo(da, NULL, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1489371c9d4SSatish Balay hx = one / (PetscReal)(mx - 1);
1499371c9d4SSatish Balay hy = one / (PetscReal)(my - 1);
1509371c9d4SSatish Balay hxdhy = hx / hy;
1519371c9d4SSatish Balay hydhx = hy / hx;
1529371c9d4SSatish Balay tleft = user->tleft;
1539371c9d4SSatish Balay tright = user->tright;
154c4762a1bSJed Brown beta = user->beta;
155c4762a1bSJed Brown
156c4762a1bSJed Brown /* Get ghost points */
1579566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
1589566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
1599566063dSJacob Faibussowitsch PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0));
1609566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArray(da, localX, &x));
1619566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArray(da, F, &f));
162c4762a1bSJed Brown
163c4762a1bSJed Brown /* Evaluate function */
164c4762a1bSJed Brown for (j = ys; j < ys + ym; j++) {
165c4762a1bSJed Brown for (i = xs; i < xs + xm; i++) {
166c4762a1bSJed Brown t0 = x[j][i];
167c4762a1bSJed Brown
168c4762a1bSJed Brown if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) {
169c4762a1bSJed Brown /* general interior volume */
170c4762a1bSJed Brown
171c4762a1bSJed Brown tw = x[j][i - 1];
172c4762a1bSJed Brown aw = 0.5 * (t0 + tw);
173c4762a1bSJed Brown dw = PetscPowScalar(aw, beta);
174c4762a1bSJed Brown fw = dw * (t0 - tw);
175c4762a1bSJed Brown
176c4762a1bSJed Brown te = x[j][i + 1];
177c4762a1bSJed Brown ae = 0.5 * (t0 + te);
178c4762a1bSJed Brown de = PetscPowScalar(ae, beta);
179c4762a1bSJed Brown fe = de * (te - t0);
180c4762a1bSJed Brown
181c4762a1bSJed Brown ts = x[j - 1][i];
182c4762a1bSJed Brown as = 0.5 * (t0 + ts);
183c4762a1bSJed Brown ds = PetscPowScalar(as, beta);
184c4762a1bSJed Brown fs = ds * (t0 - ts);
185c4762a1bSJed Brown
186c4762a1bSJed Brown tn = x[j + 1][i];
187c4762a1bSJed Brown an = 0.5 * (t0 + tn);
188c4762a1bSJed Brown dn = PetscPowScalar(an, beta);
189c4762a1bSJed Brown fn = dn * (tn - t0);
190c4762a1bSJed Brown
191c4762a1bSJed Brown } else if (i == 0) {
192c4762a1bSJed Brown /* left-hand boundary */
193c4762a1bSJed Brown tw = tleft;
194c4762a1bSJed Brown aw = 0.5 * (t0 + tw);
195c4762a1bSJed Brown dw = PetscPowScalar(aw, beta);
196c4762a1bSJed Brown fw = dw * (t0 - tw);
197c4762a1bSJed Brown
198c4762a1bSJed Brown te = x[j][i + 1];
199c4762a1bSJed Brown ae = 0.5 * (t0 + te);
200c4762a1bSJed Brown de = PetscPowScalar(ae, beta);
201c4762a1bSJed Brown fe = de * (te - t0);
202c4762a1bSJed Brown
203c4762a1bSJed Brown if (j > 0) {
204c4762a1bSJed Brown ts = x[j - 1][i];
205c4762a1bSJed Brown as = 0.5 * (t0 + ts);
206c4762a1bSJed Brown ds = PetscPowScalar(as, beta);
207c4762a1bSJed Brown fs = ds * (t0 - ts);
208c4762a1bSJed Brown } else fs = zero;
209c4762a1bSJed Brown
210c4762a1bSJed Brown if (j < my - 1) {
211c4762a1bSJed Brown tn = x[j + 1][i];
212c4762a1bSJed Brown an = 0.5 * (t0 + tn);
213c4762a1bSJed Brown dn = PetscPowScalar(an, beta);
214c4762a1bSJed Brown fn = dn * (tn - t0);
215c4762a1bSJed Brown } else fn = zero;
216c4762a1bSJed Brown
217c4762a1bSJed Brown } else if (i == mx - 1) {
218c4762a1bSJed Brown /* right-hand boundary */
219c4762a1bSJed Brown tw = x[j][i - 1];
220c4762a1bSJed Brown aw = 0.5 * (t0 + tw);
221c4762a1bSJed Brown dw = PetscPowScalar(aw, beta);
222c4762a1bSJed Brown fw = dw * (t0 - tw);
223c4762a1bSJed Brown
224c4762a1bSJed Brown te = tright;
225c4762a1bSJed Brown ae = 0.5 * (t0 + te);
226c4762a1bSJed Brown de = PetscPowScalar(ae, beta);
227c4762a1bSJed Brown fe = de * (te - t0);
228c4762a1bSJed Brown
229c4762a1bSJed Brown if (j > 0) {
230c4762a1bSJed Brown ts = x[j - 1][i];
231c4762a1bSJed Brown as = 0.5 * (t0 + ts);
232c4762a1bSJed Brown ds = PetscPowScalar(as, beta);
233c4762a1bSJed Brown fs = ds * (t0 - ts);
234c4762a1bSJed Brown } else fs = zero;
235c4762a1bSJed Brown
236c4762a1bSJed Brown if (j < my - 1) {
237c4762a1bSJed Brown tn = x[j + 1][i];
238c4762a1bSJed Brown an = 0.5 * (t0 + tn);
239c4762a1bSJed Brown dn = PetscPowScalar(an, beta);
240c4762a1bSJed Brown fn = dn * (tn - t0);
241c4762a1bSJed Brown } else fn = zero;
242c4762a1bSJed Brown
243c4762a1bSJed Brown } else if (j == 0) {
244c4762a1bSJed Brown /* bottom boundary,and i <> 0 or mx-1 */
245c4762a1bSJed Brown tw = x[j][i - 1];
246c4762a1bSJed Brown aw = 0.5 * (t0 + tw);
247c4762a1bSJed Brown dw = PetscPowScalar(aw, beta);
248c4762a1bSJed Brown fw = dw * (t0 - tw);
249c4762a1bSJed Brown
250c4762a1bSJed Brown te = x[j][i + 1];
251c4762a1bSJed Brown ae = 0.5 * (t0 + te);
252c4762a1bSJed Brown de = PetscPowScalar(ae, beta);
253c4762a1bSJed Brown fe = de * (te - t0);
254c4762a1bSJed Brown
255c4762a1bSJed Brown fs = zero;
256c4762a1bSJed Brown
257c4762a1bSJed Brown tn = x[j + 1][i];
258c4762a1bSJed Brown an = 0.5 * (t0 + tn);
259c4762a1bSJed Brown dn = PetscPowScalar(an, beta);
260c4762a1bSJed Brown fn = dn * (tn - t0);
261c4762a1bSJed Brown
262c4762a1bSJed Brown } else if (j == my - 1) {
263c4762a1bSJed Brown /* top boundary,and i <> 0 or mx-1 */
264c4762a1bSJed Brown tw = x[j][i - 1];
265c4762a1bSJed Brown aw = 0.5 * (t0 + tw);
266c4762a1bSJed Brown dw = PetscPowScalar(aw, beta);
267c4762a1bSJed Brown fw = dw * (t0 - tw);
268c4762a1bSJed Brown
269c4762a1bSJed Brown te = x[j][i + 1];
270c4762a1bSJed Brown ae = 0.5 * (t0 + te);
271c4762a1bSJed Brown de = PetscPowScalar(ae, beta);
272c4762a1bSJed Brown fe = de * (te - t0);
273c4762a1bSJed Brown
274c4762a1bSJed Brown ts = x[j - 1][i];
275c4762a1bSJed Brown as = 0.5 * (t0 + ts);
276c4762a1bSJed Brown ds = PetscPowScalar(as, beta);
277c4762a1bSJed Brown fs = ds * (t0 - ts);
278c4762a1bSJed Brown
279c4762a1bSJed Brown fn = zero;
280c4762a1bSJed Brown }
281c4762a1bSJed Brown
282c4762a1bSJed Brown f[j][i] = -hydhx * (fe - fw) - hxdhy * (fn - fs);
283c4762a1bSJed Brown }
284c4762a1bSJed Brown }
2859566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArray(da, localX, &x));
2869566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArray(da, F, &f));
2879566063dSJacob Faibussowitsch PetscCall(DMRestoreLocalVector(da, &localX));
2889566063dSJacob Faibussowitsch PetscCall(PetscLogFlops((22.0 + 4.0 * POWFLOP) * ym * xm));
2893ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
290c4762a1bSJed Brown }
291c4762a1bSJed Brown /* -------------------- Evaluate Jacobian F(x) --------------------- */
FormJacobian(SNES snes,Vec X,Mat jac,Mat B,void * ptr)292d71ae5a4SJacob Faibussowitsch PetscErrorCode FormJacobian(SNES snes, Vec X, Mat jac, Mat B, void *ptr)
293d71ae5a4SJacob Faibussowitsch {
294c4762a1bSJed Brown AppCtx *user = (AppCtx *)ptr;
295c4762a1bSJed Brown PetscInt i, j, mx, my, xs, ys, xm, ym;
296c4762a1bSJed Brown PetscScalar one = 1.0, hx, hy, hxdhy, hydhx, t0, tn, ts, te, tw;
297c4762a1bSJed Brown PetscScalar dn, ds, de, dw, an, as, ae, aw, bn, bs, be, bw, gn, gs, ge, gw;
298c4762a1bSJed Brown PetscScalar tleft, tright, beta, bm1, coef;
299c4762a1bSJed Brown PetscScalar v[5], **x;
300c4762a1bSJed Brown Vec localX;
301c4762a1bSJed Brown MatStencil col[5], row;
302c4762a1bSJed Brown DM da;
303c4762a1bSJed Brown
304c4762a1bSJed Brown PetscFunctionBeginUser;
3059566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &da));
3069566063dSJacob Faibussowitsch PetscCall(DMGetLocalVector(da, &localX));
3079566063dSJacob Faibussowitsch PetscCall(DMDAGetInfo(da, NULL, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
3089371c9d4SSatish Balay hx = one / (PetscReal)(mx - 1);
3099371c9d4SSatish Balay hy = one / (PetscReal)(my - 1);
3109371c9d4SSatish Balay hxdhy = hx / hy;
3119371c9d4SSatish Balay hydhx = hy / hx;
3129371c9d4SSatish Balay tleft = user->tleft;
3139371c9d4SSatish Balay tright = user->tright;
3149371c9d4SSatish Balay beta = user->beta;
3159371c9d4SSatish Balay bm1 = user->bm1;
3169371c9d4SSatish Balay coef = user->coef;
317c4762a1bSJed Brown
318c4762a1bSJed Brown /* Get ghost points */
3199566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
3209566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
3219566063dSJacob Faibussowitsch PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0));
3229566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArray(da, localX, &x));
323c4762a1bSJed Brown
324c4762a1bSJed Brown /* Evaluate Jacobian of function */
325c4762a1bSJed Brown for (j = ys; j < ys + ym; j++) {
326c4762a1bSJed Brown for (i = xs; i < xs + xm; i++) {
327c4762a1bSJed Brown t0 = x[j][i];
328c4762a1bSJed Brown
329c4762a1bSJed Brown if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) {
330c4762a1bSJed Brown /* general interior volume */
331c4762a1bSJed Brown
332c4762a1bSJed Brown tw = x[j][i - 1];
333c4762a1bSJed Brown aw = 0.5 * (t0 + tw);
334c4762a1bSJed Brown bw = PetscPowScalar(aw, bm1);
335c4762a1bSJed Brown /* dw = bw * aw */
336c4762a1bSJed Brown dw = PetscPowScalar(aw, beta);
337c4762a1bSJed Brown gw = coef * bw * (t0 - tw);
338c4762a1bSJed Brown
339c4762a1bSJed Brown te = x[j][i + 1];
340c4762a1bSJed Brown ae = 0.5 * (t0 + te);
341c4762a1bSJed Brown be = PetscPowScalar(ae, bm1);
342c4762a1bSJed Brown /* de = be * ae; */
343c4762a1bSJed Brown de = PetscPowScalar(ae, beta);
344c4762a1bSJed Brown ge = coef * be * (te - t0);
345c4762a1bSJed Brown
346c4762a1bSJed Brown ts = x[j - 1][i];
347c4762a1bSJed Brown as = 0.5 * (t0 + ts);
348c4762a1bSJed Brown bs = PetscPowScalar(as, bm1);
349c4762a1bSJed Brown /* ds = bs * as; */
350c4762a1bSJed Brown ds = PetscPowScalar(as, beta);
351c4762a1bSJed Brown gs = coef * bs * (t0 - ts);
352c4762a1bSJed Brown
353c4762a1bSJed Brown tn = x[j + 1][i];
354c4762a1bSJed Brown an = 0.5 * (t0 + tn);
355c4762a1bSJed Brown bn = PetscPowScalar(an, bm1);
356c4762a1bSJed Brown /* dn = bn * an; */
357c4762a1bSJed Brown dn = PetscPowScalar(an, beta);
358c4762a1bSJed Brown gn = coef * bn * (tn - t0);
359c4762a1bSJed Brown
3609371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs);
3619371c9d4SSatish Balay col[0].j = j - 1;
3629371c9d4SSatish Balay col[0].i = i;
3639371c9d4SSatish Balay v[1] = -hydhx * (dw - gw);
3649371c9d4SSatish Balay col[1].j = j;
3659371c9d4SSatish Balay col[1].i = i - 1;
3669371c9d4SSatish Balay v[2] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge);
3679371c9d4SSatish Balay col[2].j = row.j = j;
3689371c9d4SSatish Balay col[2].i = row.i = i;
3699371c9d4SSatish Balay v[3] = -hydhx * (de + ge);
3709371c9d4SSatish Balay col[3].j = j;
3719371c9d4SSatish Balay col[3].i = i + 1;
3729371c9d4SSatish Balay v[4] = -hxdhy * (dn + gn);
3739371c9d4SSatish Balay col[4].j = j + 1;
3749371c9d4SSatish Balay col[4].i = i;
3759566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 5, col, v, INSERT_VALUES));
376c4762a1bSJed Brown
377c4762a1bSJed Brown } else if (i == 0) {
378c4762a1bSJed Brown /* left-hand boundary */
379c4762a1bSJed Brown tw = tleft;
380c4762a1bSJed Brown aw = 0.5 * (t0 + tw);
381c4762a1bSJed Brown bw = PetscPowScalar(aw, bm1);
382c4762a1bSJed Brown /* dw = bw * aw */
383c4762a1bSJed Brown dw = PetscPowScalar(aw, beta);
384c4762a1bSJed Brown gw = coef * bw * (t0 - tw);
385c4762a1bSJed Brown
386c4762a1bSJed Brown te = x[j][i + 1];
387c4762a1bSJed Brown ae = 0.5 * (t0 + te);
388c4762a1bSJed Brown be = PetscPowScalar(ae, bm1);
389c4762a1bSJed Brown /* de = be * ae; */
390c4762a1bSJed Brown de = PetscPowScalar(ae, beta);
391c4762a1bSJed Brown ge = coef * be * (te - t0);
392c4762a1bSJed Brown
393c4762a1bSJed Brown /* left-hand bottom boundary */
394c4762a1bSJed Brown if (j == 0) {
395c4762a1bSJed Brown tn = x[j + 1][i];
396c4762a1bSJed Brown an = 0.5 * (t0 + tn);
397c4762a1bSJed Brown bn = PetscPowScalar(an, bm1);
398c4762a1bSJed Brown /* dn = bn * an; */
399c4762a1bSJed Brown dn = PetscPowScalar(an, beta);
400c4762a1bSJed Brown gn = coef * bn * (tn - t0);
401c4762a1bSJed Brown
4029371c9d4SSatish Balay v[0] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge);
4039371c9d4SSatish Balay col[0].j = row.j = j;
4049371c9d4SSatish Balay col[0].i = row.i = i;
4059371c9d4SSatish Balay v[1] = -hydhx * (de + ge);
4069371c9d4SSatish Balay col[1].j = j;
4079371c9d4SSatish Balay col[1].i = i + 1;
4089371c9d4SSatish Balay v[2] = -hxdhy * (dn + gn);
4099371c9d4SSatish Balay col[2].j = j + 1;
4109371c9d4SSatish Balay col[2].i = i;
4119566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES));
412c4762a1bSJed Brown
413c4762a1bSJed Brown /* left-hand interior boundary */
414c4762a1bSJed Brown } else if (j < my - 1) {
415c4762a1bSJed Brown ts = x[j - 1][i];
416c4762a1bSJed Brown as = 0.5 * (t0 + ts);
417c4762a1bSJed Brown bs = PetscPowScalar(as, bm1);
418c4762a1bSJed Brown /* ds = bs * as; */
419c4762a1bSJed Brown ds = PetscPowScalar(as, beta);
420c4762a1bSJed Brown gs = coef * bs * (t0 - ts);
421c4762a1bSJed Brown
422c4762a1bSJed Brown tn = x[j + 1][i];
423c4762a1bSJed Brown an = 0.5 * (t0 + tn);
424c4762a1bSJed Brown bn = PetscPowScalar(an, bm1);
425c4762a1bSJed Brown /* dn = bn * an; */
426c4762a1bSJed Brown dn = PetscPowScalar(an, beta);
427c4762a1bSJed Brown gn = coef * bn * (tn - t0);
428c4762a1bSJed Brown
4299371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs);
4309371c9d4SSatish Balay col[0].j = j - 1;
4319371c9d4SSatish Balay col[0].i = i;
4329371c9d4SSatish Balay v[1] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge);
4339371c9d4SSatish Balay col[1].j = row.j = j;
4349371c9d4SSatish Balay col[1].i = row.i = i;
4359371c9d4SSatish Balay v[2] = -hydhx * (de + ge);
4369371c9d4SSatish Balay col[2].j = j;
4379371c9d4SSatish Balay col[2].i = i + 1;
4389371c9d4SSatish Balay v[3] = -hxdhy * (dn + gn);
4399371c9d4SSatish Balay col[3].j = j + 1;
4409371c9d4SSatish Balay col[3].i = i;
4419566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES));
442c4762a1bSJed Brown /* left-hand top boundary */
443c4762a1bSJed Brown } else {
444c4762a1bSJed Brown ts = x[j - 1][i];
445c4762a1bSJed Brown as = 0.5 * (t0 + ts);
446c4762a1bSJed Brown bs = PetscPowScalar(as, bm1);
447c4762a1bSJed Brown /* ds = bs * as; */
448c4762a1bSJed Brown ds = PetscPowScalar(as, beta);
449c4762a1bSJed Brown gs = coef * bs * (t0 - ts);
450c4762a1bSJed Brown
4519371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs);
4529371c9d4SSatish Balay col[0].j = j - 1;
4539371c9d4SSatish Balay col[0].i = i;
4549371c9d4SSatish Balay v[1] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge);
4559371c9d4SSatish Balay col[1].j = row.j = j;
4569371c9d4SSatish Balay col[1].i = row.i = i;
4579371c9d4SSatish Balay v[2] = -hydhx * (de + ge);
4589371c9d4SSatish Balay col[2].j = j;
4599371c9d4SSatish Balay col[2].i = i + 1;
4609566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES));
461c4762a1bSJed Brown }
462c4762a1bSJed Brown
463c4762a1bSJed Brown } else if (i == mx - 1) {
464c4762a1bSJed Brown /* right-hand boundary */
465c4762a1bSJed Brown tw = x[j][i - 1];
466c4762a1bSJed Brown aw = 0.5 * (t0 + tw);
467c4762a1bSJed Brown bw = PetscPowScalar(aw, bm1);
468c4762a1bSJed Brown /* dw = bw * aw */
469c4762a1bSJed Brown dw = PetscPowScalar(aw, beta);
470c4762a1bSJed Brown gw = coef * bw * (t0 - tw);
471c4762a1bSJed Brown
472c4762a1bSJed Brown te = tright;
473c4762a1bSJed Brown ae = 0.5 * (t0 + te);
474c4762a1bSJed Brown be = PetscPowScalar(ae, bm1);
475c4762a1bSJed Brown /* de = be * ae; */
476c4762a1bSJed Brown de = PetscPowScalar(ae, beta);
477c4762a1bSJed Brown ge = coef * be * (te - t0);
478c4762a1bSJed Brown
479c4762a1bSJed Brown /* right-hand bottom boundary */
480c4762a1bSJed Brown if (j == 0) {
481c4762a1bSJed Brown tn = x[j + 1][i];
482c4762a1bSJed Brown an = 0.5 * (t0 + tn);
483c4762a1bSJed Brown bn = PetscPowScalar(an, bm1);
484c4762a1bSJed Brown /* dn = bn * an; */
485c4762a1bSJed Brown dn = PetscPowScalar(an, beta);
486c4762a1bSJed Brown gn = coef * bn * (tn - t0);
487c4762a1bSJed Brown
4889371c9d4SSatish Balay v[0] = -hydhx * (dw - gw);
4899371c9d4SSatish Balay col[0].j = j;
4909371c9d4SSatish Balay col[0].i = i - 1;
4919371c9d4SSatish Balay v[1] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge);
4929371c9d4SSatish Balay col[1].j = row.j = j;
4939371c9d4SSatish Balay col[1].i = row.i = i;
4949371c9d4SSatish Balay v[2] = -hxdhy * (dn + gn);
4959371c9d4SSatish Balay col[2].j = j + 1;
4969371c9d4SSatish Balay col[2].i = i;
4979566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES));
498c4762a1bSJed Brown
499c4762a1bSJed Brown /* right-hand interior boundary */
500c4762a1bSJed Brown } else if (j < my - 1) {
501c4762a1bSJed Brown ts = x[j - 1][i];
502c4762a1bSJed Brown as = 0.5 * (t0 + ts);
503c4762a1bSJed Brown bs = PetscPowScalar(as, bm1);
504c4762a1bSJed Brown /* ds = bs * as; */
505c4762a1bSJed Brown ds = PetscPowScalar(as, beta);
506c4762a1bSJed Brown gs = coef * bs * (t0 - ts);
507c4762a1bSJed Brown
508c4762a1bSJed Brown tn = x[j + 1][i];
509c4762a1bSJed Brown an = 0.5 * (t0 + tn);
510c4762a1bSJed Brown bn = PetscPowScalar(an, bm1);
511c4762a1bSJed Brown /* dn = bn * an; */
512c4762a1bSJed Brown dn = PetscPowScalar(an, beta);
513c4762a1bSJed Brown gn = coef * bn * (tn - t0);
514c4762a1bSJed Brown
5159371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs);
5169371c9d4SSatish Balay col[0].j = j - 1;
5179371c9d4SSatish Balay col[0].i = i;
5189371c9d4SSatish Balay v[1] = -hydhx * (dw - gw);
5199371c9d4SSatish Balay col[1].j = j;
5209371c9d4SSatish Balay col[1].i = i - 1;
5219371c9d4SSatish Balay v[2] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge);
5229371c9d4SSatish Balay col[2].j = row.j = j;
5239371c9d4SSatish Balay col[2].i = row.i = i;
5249371c9d4SSatish Balay v[3] = -hxdhy * (dn + gn);
5259371c9d4SSatish Balay col[3].j = j + 1;
5269371c9d4SSatish Balay col[3].i = i;
5279566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES));
528c4762a1bSJed Brown /* right-hand top boundary */
529c4762a1bSJed Brown } else {
530c4762a1bSJed Brown ts = x[j - 1][i];
531c4762a1bSJed Brown as = 0.5 * (t0 + ts);
532c4762a1bSJed Brown bs = PetscPowScalar(as, bm1);
533c4762a1bSJed Brown /* ds = bs * as; */
534c4762a1bSJed Brown ds = PetscPowScalar(as, beta);
535c4762a1bSJed Brown gs = coef * bs * (t0 - ts);
536c4762a1bSJed Brown
5379371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs);
5389371c9d4SSatish Balay col[0].j = j - 1;
5399371c9d4SSatish Balay col[0].i = i;
5409371c9d4SSatish Balay v[1] = -hydhx * (dw - gw);
5419371c9d4SSatish Balay col[1].j = j;
5429371c9d4SSatish Balay col[1].i = i - 1;
5439371c9d4SSatish Balay v[2] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge);
5449371c9d4SSatish Balay col[2].j = row.j = j;
5459371c9d4SSatish Balay col[2].i = row.i = i;
5469566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES));
547c4762a1bSJed Brown }
548c4762a1bSJed Brown
549c4762a1bSJed Brown /* bottom boundary,and i <> 0 or mx-1 */
550c4762a1bSJed Brown } else if (j == 0) {
551c4762a1bSJed Brown tw = x[j][i - 1];
552c4762a1bSJed Brown aw = 0.5 * (t0 + tw);
553c4762a1bSJed Brown bw = PetscPowScalar(aw, bm1);
554c4762a1bSJed Brown /* dw = bw * aw */
555c4762a1bSJed Brown dw = PetscPowScalar(aw, beta);
556c4762a1bSJed Brown gw = coef * bw * (t0 - tw);
557c4762a1bSJed Brown
558c4762a1bSJed Brown te = x[j][i + 1];
559c4762a1bSJed Brown ae = 0.5 * (t0 + te);
560c4762a1bSJed Brown be = PetscPowScalar(ae, bm1);
561c4762a1bSJed Brown /* de = be * ae; */
562c4762a1bSJed Brown de = PetscPowScalar(ae, beta);
563c4762a1bSJed Brown ge = coef * be * (te - t0);
564c4762a1bSJed Brown
565c4762a1bSJed Brown tn = x[j + 1][i];
566c4762a1bSJed Brown an = 0.5 * (t0 + tn);
567c4762a1bSJed Brown bn = PetscPowScalar(an, bm1);
568c4762a1bSJed Brown /* dn = bn * an; */
569c4762a1bSJed Brown dn = PetscPowScalar(an, beta);
570c4762a1bSJed Brown gn = coef * bn * (tn - t0);
571c4762a1bSJed Brown
5729371c9d4SSatish Balay v[0] = -hydhx * (dw - gw);
5739371c9d4SSatish Balay col[0].j = j;
5749371c9d4SSatish Balay col[0].i = i - 1;
5759371c9d4SSatish Balay v[1] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge);
5769371c9d4SSatish Balay col[1].j = row.j = j;
5779371c9d4SSatish Balay col[1].i = row.i = i;
5789371c9d4SSatish Balay v[2] = -hydhx * (de + ge);
5799371c9d4SSatish Balay col[2].j = j;
5809371c9d4SSatish Balay col[2].i = i + 1;
5819371c9d4SSatish Balay v[3] = -hxdhy * (dn + gn);
5829371c9d4SSatish Balay col[3].j = j + 1;
5839371c9d4SSatish Balay col[3].i = i;
5849566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES));
585c4762a1bSJed Brown
586c4762a1bSJed Brown /* top boundary,and i <> 0 or mx-1 */
587c4762a1bSJed Brown } else if (j == my - 1) {
588c4762a1bSJed Brown tw = x[j][i - 1];
589c4762a1bSJed Brown aw = 0.5 * (t0 + tw);
590c4762a1bSJed Brown bw = PetscPowScalar(aw, bm1);
591c4762a1bSJed Brown /* dw = bw * aw */
592c4762a1bSJed Brown dw = PetscPowScalar(aw, beta);
593c4762a1bSJed Brown gw = coef * bw * (t0 - tw);
594c4762a1bSJed Brown
595c4762a1bSJed Brown te = x[j][i + 1];
596c4762a1bSJed Brown ae = 0.5 * (t0 + te);
597c4762a1bSJed Brown be = PetscPowScalar(ae, bm1);
598c4762a1bSJed Brown /* de = be * ae; */
599c4762a1bSJed Brown de = PetscPowScalar(ae, beta);
600c4762a1bSJed Brown ge = coef * be * (te - t0);
601c4762a1bSJed Brown
602c4762a1bSJed Brown ts = x[j - 1][i];
603c4762a1bSJed Brown as = 0.5 * (t0 + ts);
604c4762a1bSJed Brown bs = PetscPowScalar(as, bm1);
605c4762a1bSJed Brown /* ds = bs * as; */
606c4762a1bSJed Brown ds = PetscPowScalar(as, beta);
607c4762a1bSJed Brown gs = coef * bs * (t0 - ts);
608c4762a1bSJed Brown
6099371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs);
6109371c9d4SSatish Balay col[0].j = j - 1;
6119371c9d4SSatish Balay col[0].i = i;
6129371c9d4SSatish Balay v[1] = -hydhx * (dw - gw);
6139371c9d4SSatish Balay col[1].j = j;
6149371c9d4SSatish Balay col[1].i = i - 1;
6159371c9d4SSatish Balay v[2] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge);
6169371c9d4SSatish Balay col[2].j = row.j = j;
6179371c9d4SSatish Balay col[2].i = row.i = i;
6189371c9d4SSatish Balay v[3] = -hydhx * (de + ge);
6199371c9d4SSatish Balay col[3].j = j;
6209371c9d4SSatish Balay col[3].i = i + 1;
6219566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES));
622c4762a1bSJed Brown }
623c4762a1bSJed Brown }
624c4762a1bSJed Brown }
6259566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
6269566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArray(da, localX, &x));
6279566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
6289566063dSJacob Faibussowitsch PetscCall(DMRestoreLocalVector(da, &localX));
629c4762a1bSJed Brown if (jac != B) {
6309566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
6319566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
632c4762a1bSJed Brown }
633c4762a1bSJed Brown
6349566063dSJacob Faibussowitsch PetscCall(PetscLogFlops((41.0 + 8.0 * POWFLOP) * xm * ym));
6353ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
636c4762a1bSJed Brown }
637c4762a1bSJed Brown
638c4762a1bSJed Brown /*TEST
639c4762a1bSJed Brown
640c4762a1bSJed Brown test:
64141ba4c6cSHeeho Park suffix: 1
642c4762a1bSJed Brown args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view
643c4762a1bSJed Brown requires: !single
644c4762a1bSJed Brown
64541ba4c6cSHeeho Park test:
64641ba4c6cSHeeho Park suffix: 2
64741ba4c6cSHeeho Park args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view -snes_type newtontrdc
64841ba4c6cSHeeho Park requires: !single
64941ba4c6cSHeeho Park
6504a221d59SStefano Zampini test:
6514a221d59SStefano Zampini suffix: 3
6524a221d59SStefano Zampini args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view -snes_type newtontr -snes_tr_fallback_type dogleg
6534a221d59SStefano Zampini requires: !single
6544a221d59SStefano Zampini
655c4762a1bSJed Brown TEST*/
656