1c4762a1bSJed Brown 2c4762a1bSJed Brown static char help[] = "Nonlinear Radiative Transport PDE with multigrid in 2d.\n\ 3c4762a1bSJed Brown Uses 2-dimensional distributed arrays.\n\ 4c4762a1bSJed Brown A 2-dim simplified Radiative Transport test problem is used, with analytic Jacobian. \n\ 5c4762a1bSJed Brown \n\ 6c4762a1bSJed Brown Solves the linear systems via multilevel methods \n\ 7c4762a1bSJed Brown \n\ 8c4762a1bSJed Brown The command line\n\ 9c4762a1bSJed Brown options are:\n\ 10c4762a1bSJed Brown -tleft <tl>, where <tl> indicates the left Diriclet BC \n\ 11c4762a1bSJed Brown -tright <tr>, where <tr> indicates the right Diriclet BC \n\ 12c4762a1bSJed Brown -beta <beta>, where <beta> indicates the exponent in T \n\n"; 13c4762a1bSJed Brown 14c4762a1bSJed Brown /* 15c4762a1bSJed Brown 16c4762a1bSJed Brown This example models the partial differential equation 17c4762a1bSJed Brown 18c4762a1bSJed Brown - Div(alpha* T^beta (GRAD T)) = 0. 19c4762a1bSJed Brown 20c4762a1bSJed Brown where beta = 2.5 and alpha = 1.0 21c4762a1bSJed Brown 22c4762a1bSJed Brown BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0. 23c4762a1bSJed Brown 24c4762a1bSJed Brown in the unit square, which is uniformly discretized in each of x and 25c4762a1bSJed Brown y in this simple encoding. The degrees of freedom are cell centered. 26c4762a1bSJed Brown 27c4762a1bSJed Brown A finite volume approximation with the usual 5-point stencil 28c4762a1bSJed Brown is used to discretize the boundary value problem to obtain a 29c4762a1bSJed Brown nonlinear system of equations. 30c4762a1bSJed Brown 31c4762a1bSJed Brown This code was contributed by David Keyes 32c4762a1bSJed Brown 33c4762a1bSJed Brown */ 34c4762a1bSJed Brown 35c4762a1bSJed Brown #include <petscsnes.h> 36c4762a1bSJed Brown #include <petscdm.h> 37c4762a1bSJed Brown #include <petscdmda.h> 38c4762a1bSJed Brown 39c4762a1bSJed Brown /* User-defined application context */ 40c4762a1bSJed Brown 41c4762a1bSJed Brown typedef struct { 42c4762a1bSJed Brown PetscReal tleft, tright; /* Dirichlet boundary conditions */ 43c4762a1bSJed Brown PetscReal beta, bm1, coef; /* nonlinear diffusivity parameterizations */ 44c4762a1bSJed Brown } AppCtx; 45c4762a1bSJed Brown 46c4762a1bSJed Brown #define POWFLOP 5 /* assume a pow() takes five flops */ 47c4762a1bSJed Brown 48c4762a1bSJed Brown extern PetscErrorCode FormInitialGuess(SNES, Vec, void *); 49c4762a1bSJed Brown extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *); 50c4762a1bSJed Brown extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *); 51c4762a1bSJed Brown 52*9371c9d4SSatish Balay int main(int argc, char **argv) { 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 ----------------- */ 108*9371c9d4SSatish Balay PetscErrorCode FormInitialGuess(SNES snes, Vec X, void *ctx) { 109c4762a1bSJed Brown AppCtx *user; 110c4762a1bSJed Brown PetscInt i, j, xs, ys, xm, ym; 111c4762a1bSJed Brown PetscReal tleft; 112c4762a1bSJed Brown PetscScalar **x; 113c4762a1bSJed Brown DM da; 114c4762a1bSJed Brown 115c4762a1bSJed Brown PetscFunctionBeginUser; 1169566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &da)); 1179566063dSJacob Faibussowitsch PetscCall(DMGetApplicationContext(da, &user)); 118c4762a1bSJed Brown tleft = user->tleft; 119c4762a1bSJed Brown /* Get ghost points */ 1209566063dSJacob Faibussowitsch PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0)); 1219566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArray(da, X, &x)); 122c4762a1bSJed Brown 123c4762a1bSJed Brown /* Compute initial guess */ 124c4762a1bSJed Brown for (j = ys; j < ys + ym; j++) { 125*9371c9d4SSatish Balay for (i = xs; i < xs + xm; i++) { x[j][i] = tleft; } 126c4762a1bSJed Brown } 1279566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArray(da, X, &x)); 128c4762a1bSJed Brown PetscFunctionReturn(0); 129c4762a1bSJed Brown } 130c4762a1bSJed Brown /* -------------------- Evaluate Function F(x) --------------------- */ 131*9371c9d4SSatish Balay PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr) { 132c4762a1bSJed Brown AppCtx *user = (AppCtx *)ptr; 133c4762a1bSJed Brown PetscInt i, j, mx, my, xs, ys, xm, ym; 134c4762a1bSJed Brown PetscScalar zero = 0.0, one = 1.0; 135c4762a1bSJed Brown PetscScalar hx, hy, hxdhy, hydhx; 136c4762a1bSJed 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; 137c4762a1bSJed Brown PetscScalar tleft, tright, beta; 138c4762a1bSJed Brown PetscScalar **x, **f; 139c4762a1bSJed Brown Vec localX; 140c4762a1bSJed Brown DM da; 141c4762a1bSJed Brown 142c4762a1bSJed Brown PetscFunctionBeginUser; 1439566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &da)); 1449566063dSJacob Faibussowitsch PetscCall(DMGetLocalVector(da, &localX)); 1459566063dSJacob Faibussowitsch PetscCall(DMDAGetInfo(da, NULL, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 146*9371c9d4SSatish Balay hx = one / (PetscReal)(mx - 1); 147*9371c9d4SSatish Balay hy = one / (PetscReal)(my - 1); 148*9371c9d4SSatish Balay hxdhy = hx / hy; 149*9371c9d4SSatish Balay hydhx = hy / hx; 150*9371c9d4SSatish Balay tleft = user->tleft; 151*9371c9d4SSatish Balay tright = user->tright; 152c4762a1bSJed Brown beta = user->beta; 153c4762a1bSJed Brown 154c4762a1bSJed Brown /* Get ghost points */ 1559566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 1569566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 1579566063dSJacob Faibussowitsch PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0)); 1589566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArray(da, localX, &x)); 1599566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArray(da, F, &f)); 160c4762a1bSJed Brown 161c4762a1bSJed Brown /* Evaluate function */ 162c4762a1bSJed Brown for (j = ys; j < ys + ym; j++) { 163c4762a1bSJed Brown for (i = xs; i < xs + xm; i++) { 164c4762a1bSJed Brown t0 = x[j][i]; 165c4762a1bSJed Brown 166c4762a1bSJed Brown if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) { 167c4762a1bSJed Brown /* general interior volume */ 168c4762a1bSJed Brown 169c4762a1bSJed Brown tw = x[j][i - 1]; 170c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 171c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 172c4762a1bSJed Brown fw = dw * (t0 - tw); 173c4762a1bSJed Brown 174c4762a1bSJed Brown te = x[j][i + 1]; 175c4762a1bSJed Brown ae = 0.5 * (t0 + te); 176c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 177c4762a1bSJed Brown fe = de * (te - t0); 178c4762a1bSJed Brown 179c4762a1bSJed Brown ts = x[j - 1][i]; 180c4762a1bSJed Brown as = 0.5 * (t0 + ts); 181c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 182c4762a1bSJed Brown fs = ds * (t0 - ts); 183c4762a1bSJed Brown 184c4762a1bSJed Brown tn = x[j + 1][i]; 185c4762a1bSJed Brown an = 0.5 * (t0 + tn); 186c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 187c4762a1bSJed Brown fn = dn * (tn - t0); 188c4762a1bSJed Brown 189c4762a1bSJed Brown } else if (i == 0) { 190c4762a1bSJed Brown /* left-hand boundary */ 191c4762a1bSJed Brown tw = tleft; 192c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 193c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 194c4762a1bSJed Brown fw = dw * (t0 - tw); 195c4762a1bSJed Brown 196c4762a1bSJed Brown te = x[j][i + 1]; 197c4762a1bSJed Brown ae = 0.5 * (t0 + te); 198c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 199c4762a1bSJed Brown fe = de * (te - t0); 200c4762a1bSJed Brown 201c4762a1bSJed Brown if (j > 0) { 202c4762a1bSJed Brown ts = x[j - 1][i]; 203c4762a1bSJed Brown as = 0.5 * (t0 + ts); 204c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 205c4762a1bSJed Brown fs = ds * (t0 - ts); 206c4762a1bSJed Brown } else fs = zero; 207c4762a1bSJed Brown 208c4762a1bSJed Brown if (j < my - 1) { 209c4762a1bSJed Brown tn = x[j + 1][i]; 210c4762a1bSJed Brown an = 0.5 * (t0 + tn); 211c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 212c4762a1bSJed Brown fn = dn * (tn - t0); 213c4762a1bSJed Brown } else fn = zero; 214c4762a1bSJed Brown 215c4762a1bSJed Brown } else if (i == mx - 1) { 216c4762a1bSJed Brown /* right-hand boundary */ 217c4762a1bSJed Brown tw = x[j][i - 1]; 218c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 219c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 220c4762a1bSJed Brown fw = dw * (t0 - tw); 221c4762a1bSJed Brown 222c4762a1bSJed Brown te = tright; 223c4762a1bSJed Brown ae = 0.5 * (t0 + te); 224c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 225c4762a1bSJed Brown fe = de * (te - t0); 226c4762a1bSJed Brown 227c4762a1bSJed Brown if (j > 0) { 228c4762a1bSJed Brown ts = x[j - 1][i]; 229c4762a1bSJed Brown as = 0.5 * (t0 + ts); 230c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 231c4762a1bSJed Brown fs = ds * (t0 - ts); 232c4762a1bSJed Brown } else fs = zero; 233c4762a1bSJed Brown 234c4762a1bSJed Brown if (j < my - 1) { 235c4762a1bSJed Brown tn = x[j + 1][i]; 236c4762a1bSJed Brown an = 0.5 * (t0 + tn); 237c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 238c4762a1bSJed Brown fn = dn * (tn - t0); 239c4762a1bSJed Brown } else fn = zero; 240c4762a1bSJed Brown 241c4762a1bSJed Brown } else if (j == 0) { 242c4762a1bSJed Brown /* bottom boundary,and i <> 0 or mx-1 */ 243c4762a1bSJed Brown tw = x[j][i - 1]; 244c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 245c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 246c4762a1bSJed Brown fw = dw * (t0 - tw); 247c4762a1bSJed Brown 248c4762a1bSJed Brown te = x[j][i + 1]; 249c4762a1bSJed Brown ae = 0.5 * (t0 + te); 250c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 251c4762a1bSJed Brown fe = de * (te - t0); 252c4762a1bSJed Brown 253c4762a1bSJed Brown fs = zero; 254c4762a1bSJed Brown 255c4762a1bSJed Brown tn = x[j + 1][i]; 256c4762a1bSJed Brown an = 0.5 * (t0 + tn); 257c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 258c4762a1bSJed Brown fn = dn * (tn - t0); 259c4762a1bSJed Brown 260c4762a1bSJed Brown } else if (j == my - 1) { 261c4762a1bSJed Brown /* top boundary,and i <> 0 or mx-1 */ 262c4762a1bSJed Brown tw = x[j][i - 1]; 263c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 264c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 265c4762a1bSJed Brown fw = dw * (t0 - tw); 266c4762a1bSJed Brown 267c4762a1bSJed Brown te = x[j][i + 1]; 268c4762a1bSJed Brown ae = 0.5 * (t0 + te); 269c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 270c4762a1bSJed Brown fe = de * (te - t0); 271c4762a1bSJed Brown 272c4762a1bSJed Brown ts = x[j - 1][i]; 273c4762a1bSJed Brown as = 0.5 * (t0 + ts); 274c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 275c4762a1bSJed Brown fs = ds * (t0 - ts); 276c4762a1bSJed Brown 277c4762a1bSJed Brown fn = zero; 278c4762a1bSJed Brown } 279c4762a1bSJed Brown 280c4762a1bSJed Brown f[j][i] = -hydhx * (fe - fw) - hxdhy * (fn - fs); 281c4762a1bSJed Brown } 282c4762a1bSJed Brown } 2839566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArray(da, localX, &x)); 2849566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArray(da, F, &f)); 2859566063dSJacob Faibussowitsch PetscCall(DMRestoreLocalVector(da, &localX)); 2869566063dSJacob Faibussowitsch PetscCall(PetscLogFlops((22.0 + 4.0 * POWFLOP) * ym * xm)); 287c4762a1bSJed Brown PetscFunctionReturn(0); 288c4762a1bSJed Brown } 289c4762a1bSJed Brown /* -------------------- Evaluate Jacobian F(x) --------------------- */ 290*9371c9d4SSatish Balay PetscErrorCode FormJacobian(SNES snes, Vec X, Mat jac, Mat B, void *ptr) { 291c4762a1bSJed Brown AppCtx *user = (AppCtx *)ptr; 292c4762a1bSJed Brown PetscInt i, j, mx, my, xs, ys, xm, ym; 293c4762a1bSJed Brown PetscScalar one = 1.0, hx, hy, hxdhy, hydhx, t0, tn, ts, te, tw; 294c4762a1bSJed Brown PetscScalar dn, ds, de, dw, an, as, ae, aw, bn, bs, be, bw, gn, gs, ge, gw; 295c4762a1bSJed Brown PetscScalar tleft, tright, beta, bm1, coef; 296c4762a1bSJed Brown PetscScalar v[5], **x; 297c4762a1bSJed Brown Vec localX; 298c4762a1bSJed Brown MatStencil col[5], row; 299c4762a1bSJed Brown DM da; 300c4762a1bSJed Brown 301c4762a1bSJed Brown PetscFunctionBeginUser; 3029566063dSJacob Faibussowitsch PetscCall(SNESGetDM(snes, &da)); 3039566063dSJacob Faibussowitsch PetscCall(DMGetLocalVector(da, &localX)); 3049566063dSJacob Faibussowitsch PetscCall(DMDAGetInfo(da, NULL, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 305*9371c9d4SSatish Balay hx = one / (PetscReal)(mx - 1); 306*9371c9d4SSatish Balay hy = one / (PetscReal)(my - 1); 307*9371c9d4SSatish Balay hxdhy = hx / hy; 308*9371c9d4SSatish Balay hydhx = hy / hx; 309*9371c9d4SSatish Balay tleft = user->tleft; 310*9371c9d4SSatish Balay tright = user->tright; 311*9371c9d4SSatish Balay beta = user->beta; 312*9371c9d4SSatish Balay bm1 = user->bm1; 313*9371c9d4SSatish Balay coef = user->coef; 314c4762a1bSJed Brown 315c4762a1bSJed Brown /* Get ghost points */ 3169566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 3179566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 3189566063dSJacob Faibussowitsch PetscCall(DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0)); 3199566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArray(da, localX, &x)); 320c4762a1bSJed Brown 321c4762a1bSJed Brown /* Evaluate Jacobian of function */ 322c4762a1bSJed Brown for (j = ys; j < ys + ym; j++) { 323c4762a1bSJed Brown for (i = xs; i < xs + xm; i++) { 324c4762a1bSJed Brown t0 = x[j][i]; 325c4762a1bSJed Brown 326c4762a1bSJed Brown if (i > 0 && i < mx - 1 && j > 0 && j < my - 1) { 327c4762a1bSJed Brown /* general interior volume */ 328c4762a1bSJed Brown 329c4762a1bSJed Brown tw = x[j][i - 1]; 330c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 331c4762a1bSJed Brown bw = PetscPowScalar(aw, bm1); 332c4762a1bSJed Brown /* dw = bw * aw */ 333c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 334c4762a1bSJed Brown gw = coef * bw * (t0 - tw); 335c4762a1bSJed Brown 336c4762a1bSJed Brown te = x[j][i + 1]; 337c4762a1bSJed Brown ae = 0.5 * (t0 + te); 338c4762a1bSJed Brown be = PetscPowScalar(ae, bm1); 339c4762a1bSJed Brown /* de = be * ae; */ 340c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 341c4762a1bSJed Brown ge = coef * be * (te - t0); 342c4762a1bSJed Brown 343c4762a1bSJed Brown ts = x[j - 1][i]; 344c4762a1bSJed Brown as = 0.5 * (t0 + ts); 345c4762a1bSJed Brown bs = PetscPowScalar(as, bm1); 346c4762a1bSJed Brown /* ds = bs * as; */ 347c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 348c4762a1bSJed Brown gs = coef * bs * (t0 - ts); 349c4762a1bSJed Brown 350c4762a1bSJed Brown tn = x[j + 1][i]; 351c4762a1bSJed Brown an = 0.5 * (t0 + tn); 352c4762a1bSJed Brown bn = PetscPowScalar(an, bm1); 353c4762a1bSJed Brown /* dn = bn * an; */ 354c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 355c4762a1bSJed Brown gn = coef * bn * (tn - t0); 356c4762a1bSJed Brown 357*9371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs); 358*9371c9d4SSatish Balay col[0].j = j - 1; 359*9371c9d4SSatish Balay col[0].i = i; 360*9371c9d4SSatish Balay v[1] = -hydhx * (dw - gw); 361*9371c9d4SSatish Balay col[1].j = j; 362*9371c9d4SSatish Balay col[1].i = i - 1; 363*9371c9d4SSatish Balay v[2] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge); 364*9371c9d4SSatish Balay col[2].j = row.j = j; 365*9371c9d4SSatish Balay col[2].i = row.i = i; 366*9371c9d4SSatish Balay v[3] = -hydhx * (de + ge); 367*9371c9d4SSatish Balay col[3].j = j; 368*9371c9d4SSatish Balay col[3].i = i + 1; 369*9371c9d4SSatish Balay v[4] = -hxdhy * (dn + gn); 370*9371c9d4SSatish Balay col[4].j = j + 1; 371*9371c9d4SSatish Balay col[4].i = i; 3729566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 5, col, v, INSERT_VALUES)); 373c4762a1bSJed Brown 374c4762a1bSJed Brown } else if (i == 0) { 375c4762a1bSJed Brown /* left-hand boundary */ 376c4762a1bSJed Brown tw = tleft; 377c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 378c4762a1bSJed Brown bw = PetscPowScalar(aw, bm1); 379c4762a1bSJed Brown /* dw = bw * aw */ 380c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 381c4762a1bSJed Brown gw = coef * bw * (t0 - tw); 382c4762a1bSJed Brown 383c4762a1bSJed Brown te = x[j][i + 1]; 384c4762a1bSJed Brown ae = 0.5 * (t0 + te); 385c4762a1bSJed Brown be = PetscPowScalar(ae, bm1); 386c4762a1bSJed Brown /* de = be * ae; */ 387c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 388c4762a1bSJed Brown ge = coef * be * (te - t0); 389c4762a1bSJed Brown 390c4762a1bSJed Brown /* left-hand bottom boundary */ 391c4762a1bSJed Brown if (j == 0) { 392c4762a1bSJed Brown tn = x[j + 1][i]; 393c4762a1bSJed Brown an = 0.5 * (t0 + tn); 394c4762a1bSJed Brown bn = PetscPowScalar(an, bm1); 395c4762a1bSJed Brown /* dn = bn * an; */ 396c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 397c4762a1bSJed Brown gn = coef * bn * (tn - t0); 398c4762a1bSJed Brown 399*9371c9d4SSatish Balay v[0] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge); 400*9371c9d4SSatish Balay col[0].j = row.j = j; 401*9371c9d4SSatish Balay col[0].i = row.i = i; 402*9371c9d4SSatish Balay v[1] = -hydhx * (de + ge); 403*9371c9d4SSatish Balay col[1].j = j; 404*9371c9d4SSatish Balay col[1].i = i + 1; 405*9371c9d4SSatish Balay v[2] = -hxdhy * (dn + gn); 406*9371c9d4SSatish Balay col[2].j = j + 1; 407*9371c9d4SSatish Balay col[2].i = i; 4089566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES)); 409c4762a1bSJed Brown 410c4762a1bSJed Brown /* left-hand interior boundary */ 411c4762a1bSJed Brown } else if (j < my - 1) { 412c4762a1bSJed Brown ts = x[j - 1][i]; 413c4762a1bSJed Brown as = 0.5 * (t0 + ts); 414c4762a1bSJed Brown bs = PetscPowScalar(as, bm1); 415c4762a1bSJed Brown /* ds = bs * as; */ 416c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 417c4762a1bSJed Brown gs = coef * bs * (t0 - ts); 418c4762a1bSJed Brown 419c4762a1bSJed Brown tn = x[j + 1][i]; 420c4762a1bSJed Brown an = 0.5 * (t0 + tn); 421c4762a1bSJed Brown bn = PetscPowScalar(an, bm1); 422c4762a1bSJed Brown /* dn = bn * an; */ 423c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 424c4762a1bSJed Brown gn = coef * bn * (tn - t0); 425c4762a1bSJed Brown 426*9371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs); 427*9371c9d4SSatish Balay col[0].j = j - 1; 428*9371c9d4SSatish Balay col[0].i = i; 429*9371c9d4SSatish Balay v[1] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge); 430*9371c9d4SSatish Balay col[1].j = row.j = j; 431*9371c9d4SSatish Balay col[1].i = row.i = i; 432*9371c9d4SSatish Balay v[2] = -hydhx * (de + ge); 433*9371c9d4SSatish Balay col[2].j = j; 434*9371c9d4SSatish Balay col[2].i = i + 1; 435*9371c9d4SSatish Balay v[3] = -hxdhy * (dn + gn); 436*9371c9d4SSatish Balay col[3].j = j + 1; 437*9371c9d4SSatish Balay col[3].i = i; 4389566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES)); 439c4762a1bSJed Brown /* left-hand top boundary */ 440c4762a1bSJed Brown } else { 441c4762a1bSJed Brown ts = x[j - 1][i]; 442c4762a1bSJed Brown as = 0.5 * (t0 + ts); 443c4762a1bSJed Brown bs = PetscPowScalar(as, bm1); 444c4762a1bSJed Brown /* ds = bs * as; */ 445c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 446c4762a1bSJed Brown gs = coef * bs * (t0 - ts); 447c4762a1bSJed Brown 448*9371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs); 449*9371c9d4SSatish Balay col[0].j = j - 1; 450*9371c9d4SSatish Balay col[0].i = i; 451*9371c9d4SSatish Balay v[1] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge); 452*9371c9d4SSatish Balay col[1].j = row.j = j; 453*9371c9d4SSatish Balay col[1].i = row.i = i; 454*9371c9d4SSatish Balay v[2] = -hydhx * (de + ge); 455*9371c9d4SSatish Balay col[2].j = j; 456*9371c9d4SSatish Balay col[2].i = i + 1; 4579566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES)); 458c4762a1bSJed Brown } 459c4762a1bSJed Brown 460c4762a1bSJed Brown } else if (i == mx - 1) { 461c4762a1bSJed Brown /* right-hand boundary */ 462c4762a1bSJed Brown tw = x[j][i - 1]; 463c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 464c4762a1bSJed Brown bw = PetscPowScalar(aw, bm1); 465c4762a1bSJed Brown /* dw = bw * aw */ 466c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 467c4762a1bSJed Brown gw = coef * bw * (t0 - tw); 468c4762a1bSJed Brown 469c4762a1bSJed Brown te = tright; 470c4762a1bSJed Brown ae = 0.5 * (t0 + te); 471c4762a1bSJed Brown be = PetscPowScalar(ae, bm1); 472c4762a1bSJed Brown /* de = be * ae; */ 473c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 474c4762a1bSJed Brown ge = coef * be * (te - t0); 475c4762a1bSJed Brown 476c4762a1bSJed Brown /* right-hand bottom boundary */ 477c4762a1bSJed Brown if (j == 0) { 478c4762a1bSJed Brown tn = x[j + 1][i]; 479c4762a1bSJed Brown an = 0.5 * (t0 + tn); 480c4762a1bSJed Brown bn = PetscPowScalar(an, bm1); 481c4762a1bSJed Brown /* dn = bn * an; */ 482c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 483c4762a1bSJed Brown gn = coef * bn * (tn - t0); 484c4762a1bSJed Brown 485*9371c9d4SSatish Balay v[0] = -hydhx * (dw - gw); 486*9371c9d4SSatish Balay col[0].j = j; 487*9371c9d4SSatish Balay col[0].i = i - 1; 488*9371c9d4SSatish Balay v[1] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge); 489*9371c9d4SSatish Balay col[1].j = row.j = j; 490*9371c9d4SSatish Balay col[1].i = row.i = i; 491*9371c9d4SSatish Balay v[2] = -hxdhy * (dn + gn); 492*9371c9d4SSatish Balay col[2].j = j + 1; 493*9371c9d4SSatish Balay col[2].i = i; 4949566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES)); 495c4762a1bSJed Brown 496c4762a1bSJed Brown /* right-hand interior boundary */ 497c4762a1bSJed Brown } else if (j < my - 1) { 498c4762a1bSJed Brown ts = x[j - 1][i]; 499c4762a1bSJed Brown as = 0.5 * (t0 + ts); 500c4762a1bSJed Brown bs = PetscPowScalar(as, bm1); 501c4762a1bSJed Brown /* ds = bs * as; */ 502c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 503c4762a1bSJed Brown gs = coef * bs * (t0 - ts); 504c4762a1bSJed Brown 505c4762a1bSJed Brown tn = x[j + 1][i]; 506c4762a1bSJed Brown an = 0.5 * (t0 + tn); 507c4762a1bSJed Brown bn = PetscPowScalar(an, bm1); 508c4762a1bSJed Brown /* dn = bn * an; */ 509c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 510c4762a1bSJed Brown gn = coef * bn * (tn - t0); 511c4762a1bSJed Brown 512*9371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs); 513*9371c9d4SSatish Balay col[0].j = j - 1; 514*9371c9d4SSatish Balay col[0].i = i; 515*9371c9d4SSatish Balay v[1] = -hydhx * (dw - gw); 516*9371c9d4SSatish Balay col[1].j = j; 517*9371c9d4SSatish Balay col[1].i = i - 1; 518*9371c9d4SSatish Balay v[2] = hxdhy * (ds + dn + gs - gn) + hydhx * (dw + de + gw - ge); 519*9371c9d4SSatish Balay col[2].j = row.j = j; 520*9371c9d4SSatish Balay col[2].i = row.i = i; 521*9371c9d4SSatish Balay v[3] = -hxdhy * (dn + gn); 522*9371c9d4SSatish Balay col[3].j = j + 1; 523*9371c9d4SSatish Balay col[3].i = i; 5249566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES)); 525c4762a1bSJed Brown /* right-hand top boundary */ 526c4762a1bSJed Brown } else { 527c4762a1bSJed Brown ts = x[j - 1][i]; 528c4762a1bSJed Brown as = 0.5 * (t0 + ts); 529c4762a1bSJed Brown bs = PetscPowScalar(as, bm1); 530c4762a1bSJed Brown /* ds = bs * as; */ 531c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 532c4762a1bSJed Brown gs = coef * bs * (t0 - ts); 533c4762a1bSJed Brown 534*9371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs); 535*9371c9d4SSatish Balay col[0].j = j - 1; 536*9371c9d4SSatish Balay col[0].i = i; 537*9371c9d4SSatish Balay v[1] = -hydhx * (dw - gw); 538*9371c9d4SSatish Balay col[1].j = j; 539*9371c9d4SSatish Balay col[1].i = i - 1; 540*9371c9d4SSatish Balay v[2] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge); 541*9371c9d4SSatish Balay col[2].j = row.j = j; 542*9371c9d4SSatish Balay col[2].i = row.i = i; 5439566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 3, col, v, INSERT_VALUES)); 544c4762a1bSJed Brown } 545c4762a1bSJed Brown 546c4762a1bSJed Brown /* bottom boundary,and i <> 0 or mx-1 */ 547c4762a1bSJed Brown } else if (j == 0) { 548c4762a1bSJed Brown tw = x[j][i - 1]; 549c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 550c4762a1bSJed Brown bw = PetscPowScalar(aw, bm1); 551c4762a1bSJed Brown /* dw = bw * aw */ 552c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 553c4762a1bSJed Brown gw = coef * bw * (t0 - tw); 554c4762a1bSJed Brown 555c4762a1bSJed Brown te = x[j][i + 1]; 556c4762a1bSJed Brown ae = 0.5 * (t0 + te); 557c4762a1bSJed Brown be = PetscPowScalar(ae, bm1); 558c4762a1bSJed Brown /* de = be * ae; */ 559c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 560c4762a1bSJed Brown ge = coef * be * (te - t0); 561c4762a1bSJed Brown 562c4762a1bSJed Brown tn = x[j + 1][i]; 563c4762a1bSJed Brown an = 0.5 * (t0 + tn); 564c4762a1bSJed Brown bn = PetscPowScalar(an, bm1); 565c4762a1bSJed Brown /* dn = bn * an; */ 566c4762a1bSJed Brown dn = PetscPowScalar(an, beta); 567c4762a1bSJed Brown gn = coef * bn * (tn - t0); 568c4762a1bSJed Brown 569*9371c9d4SSatish Balay v[0] = -hydhx * (dw - gw); 570*9371c9d4SSatish Balay col[0].j = j; 571*9371c9d4SSatish Balay col[0].i = i - 1; 572*9371c9d4SSatish Balay v[1] = hxdhy * (dn - gn) + hydhx * (dw + de + gw - ge); 573*9371c9d4SSatish Balay col[1].j = row.j = j; 574*9371c9d4SSatish Balay col[1].i = row.i = i; 575*9371c9d4SSatish Balay v[2] = -hydhx * (de + ge); 576*9371c9d4SSatish Balay col[2].j = j; 577*9371c9d4SSatish Balay col[2].i = i + 1; 578*9371c9d4SSatish Balay v[3] = -hxdhy * (dn + gn); 579*9371c9d4SSatish Balay col[3].j = j + 1; 580*9371c9d4SSatish Balay col[3].i = i; 5819566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES)); 582c4762a1bSJed Brown 583c4762a1bSJed Brown /* top boundary,and i <> 0 or mx-1 */ 584c4762a1bSJed Brown } else if (j == my - 1) { 585c4762a1bSJed Brown tw = x[j][i - 1]; 586c4762a1bSJed Brown aw = 0.5 * (t0 + tw); 587c4762a1bSJed Brown bw = PetscPowScalar(aw, bm1); 588c4762a1bSJed Brown /* dw = bw * aw */ 589c4762a1bSJed Brown dw = PetscPowScalar(aw, beta); 590c4762a1bSJed Brown gw = coef * bw * (t0 - tw); 591c4762a1bSJed Brown 592c4762a1bSJed Brown te = x[j][i + 1]; 593c4762a1bSJed Brown ae = 0.5 * (t0 + te); 594c4762a1bSJed Brown be = PetscPowScalar(ae, bm1); 595c4762a1bSJed Brown /* de = be * ae; */ 596c4762a1bSJed Brown de = PetscPowScalar(ae, beta); 597c4762a1bSJed Brown ge = coef * be * (te - t0); 598c4762a1bSJed Brown 599c4762a1bSJed Brown ts = x[j - 1][i]; 600c4762a1bSJed Brown as = 0.5 * (t0 + ts); 601c4762a1bSJed Brown bs = PetscPowScalar(as, bm1); 602c4762a1bSJed Brown /* ds = bs * as; */ 603c4762a1bSJed Brown ds = PetscPowScalar(as, beta); 604c4762a1bSJed Brown gs = coef * bs * (t0 - ts); 605c4762a1bSJed Brown 606*9371c9d4SSatish Balay v[0] = -hxdhy * (ds - gs); 607*9371c9d4SSatish Balay col[0].j = j - 1; 608*9371c9d4SSatish Balay col[0].i = i; 609*9371c9d4SSatish Balay v[1] = -hydhx * (dw - gw); 610*9371c9d4SSatish Balay col[1].j = j; 611*9371c9d4SSatish Balay col[1].i = i - 1; 612*9371c9d4SSatish Balay v[2] = hxdhy * (ds + gs) + hydhx * (dw + de + gw - ge); 613*9371c9d4SSatish Balay col[2].j = row.j = j; 614*9371c9d4SSatish Balay col[2].i = row.i = i; 615*9371c9d4SSatish Balay v[3] = -hydhx * (de + ge); 616*9371c9d4SSatish Balay col[3].j = j; 617*9371c9d4SSatish Balay col[3].i = i + 1; 6189566063dSJacob Faibussowitsch PetscCall(MatSetValuesStencil(B, 1, &row, 4, col, v, INSERT_VALUES)); 619c4762a1bSJed Brown } 620c4762a1bSJed Brown } 621c4762a1bSJed Brown } 6229566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 6239566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArray(da, localX, &x)); 6249566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 6259566063dSJacob Faibussowitsch PetscCall(DMRestoreLocalVector(da, &localX)); 626c4762a1bSJed Brown if (jac != B) { 6279566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 6289566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 629c4762a1bSJed Brown } 630c4762a1bSJed Brown 6319566063dSJacob Faibussowitsch PetscCall(PetscLogFlops((41.0 + 8.0 * POWFLOP) * xm * ym)); 632c4762a1bSJed Brown PetscFunctionReturn(0); 633c4762a1bSJed Brown } 634c4762a1bSJed Brown 635c4762a1bSJed Brown /*TEST 636c4762a1bSJed Brown 637c4762a1bSJed Brown test: 63841ba4c6cSHeeho Park suffix: 1 639c4762a1bSJed Brown args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view 640c4762a1bSJed Brown requires: !single 641c4762a1bSJed Brown 64241ba4c6cSHeeho Park test: 64341ba4c6cSHeeho Park suffix: 2 64441ba4c6cSHeeho Park args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view -snes_type newtontrdc 64541ba4c6cSHeeho Park requires: !single 64641ba4c6cSHeeho Park 647c4762a1bSJed Brown TEST*/ 648