Actual source code: ex10.c

  1: /*
  2:   This example implements the model described in

  4:     Rauenzahn, Mousseau, Knoll. "Temporal accuracy of the nonequilibrium radiation diffusion
  5:     equations employing a Saha ionization model" 2005.

  7:   The paper discusses three examples, the first two are nondimensional with a simple
  8:   ionization model.  The third example is fully dimensional and uses the Saha ionization
  9:   model with realistic parameters.
 10: */

 12: #include <petscts.h>
 13: #include <petscdm.h>
 14: #include <petscdmda.h>

 16: typedef enum {
 17:   BC_DIRICHLET,
 18:   BC_NEUMANN,
 19:   BC_ROBIN
 20: } BCType;
 21: static const char *const BCTypes[] = {"DIRICHLET", "NEUMANN", "ROBIN", "BCType", "BC_", 0};
 22: typedef enum {
 23:   JACOBIAN_ANALYTIC,
 24:   JACOBIAN_MATRIXFREE,
 25:   JACOBIAN_FD_COLORING,
 26:   JACOBIAN_FD_FULL
 27: } JacobianType;
 28: static const char *const JacobianTypes[] = {"ANALYTIC", "MATRIXFREE", "FD_COLORING", "FD_FULL", "JacobianType", "FD_", 0};
 29: typedef enum {
 30:   DISCRETIZATION_FD,
 31:   DISCRETIZATION_FE
 32: } DiscretizationType;
 33: static const char *const DiscretizationTypes[] = {"FD", "FE", "DiscretizationType", "DISCRETIZATION_", 0};
 34: typedef enum {
 35:   QUADRATURE_GAUSS1,
 36:   QUADRATURE_GAUSS2,
 37:   QUADRATURE_GAUSS3,
 38:   QUADRATURE_GAUSS4,
 39:   QUADRATURE_LOBATTO2,
 40:   QUADRATURE_LOBATTO3
 41: } QuadratureType;
 42: static const char *const QuadratureTypes[] = {"GAUSS1", "GAUSS2", "GAUSS3", "GAUSS4", "LOBATTO2", "LOBATTO3", "QuadratureType", "QUADRATURE_", 0};

 44: typedef struct {
 45:   PetscScalar E; /* radiation energy */
 46:   PetscScalar T; /* material temperature */
 47: } RDNode;

 49: typedef struct {
 50:   PetscReal meter, kilogram, second, Kelvin; /* Fundamental units */
 51:   PetscReal Joule, Watt;                     /* Derived units */
 52: } RDUnit;

 54: typedef struct _n_RD *RD;

 56: struct _n_RD {
 57:   void (*MaterialEnergy)(RD, const RDNode *, PetscScalar *, RDNode *);
 58:   DM                 da;
 59:   PetscBool          monitor_residual;
 60:   DiscretizationType discretization;
 61:   QuadratureType     quadrature;
 62:   JacobianType       jacobian;
 63:   PetscInt           initial;
 64:   BCType             leftbc;
 65:   PetscBool          view_draw;
 66:   char               view_binary[PETSC_MAX_PATH_LEN];
 67:   PetscBool          test_diff;
 68:   PetscBool          endpoint;
 69:   PetscBool          bclimit;
 70:   PetscBool          bcmidpoint;
 71:   RDUnit             unit;

 73:   /* model constants, see Table 2 and RDCreate() */
 74:   PetscReal rho, K_R, K_p, I_H, m_p, m_e, h, k, c, sigma_b, beta, gamma;

 76:   /* Domain and boundary conditions */
 77:   PetscReal Eapplied; /* Radiation flux from the left */
 78:   PetscReal L;        /* Length of domain */
 79:   PetscReal final_time;
 80: };

 82: static PetscErrorCode RDDestroy(RD *rd)
 83: {
 85:   DMDestroy(&(*rd)->da);
 86:   PetscFree(*rd);
 87:   return 0;
 88: }

 90: /* The paper has a time derivative for material energy (Eq 2) which is a dependent variable (computable from temperature
 91:  * and density through an uninvertible relation).  Computing this derivative is trivial for trapezoid rule (used in the
 92:  * paper), but does not generalize nicely to higher order integrators.  Here we use the implicit form which provides
 93:  * time derivatives of the independent variables (radiation energy and temperature), so we must compute the time
 94:  * derivative of material energy ourselves (could be done using AD).
 95:  *
 96:  * There are multiple ionization models, this interface dispatches to the one currently in use.
 97:  */
 98: static void RDMaterialEnergy(RD rd, const RDNode *n, PetscScalar *Em, RDNode *dEm)
 99: {
100:   rd->MaterialEnergy(rd, n, Em, dEm);
101: }

103: /* Solves a quadratic equation while propagating tangents */
104: static void QuadraticSolve(PetscScalar a, PetscScalar a_t, PetscScalar b, PetscScalar b_t, PetscScalar c, PetscScalar c_t, PetscScalar *x, PetscScalar *x_t)
105: {
106:   PetscScalar disc = b * b - 4. * a * c, disc_t = 2. * b * b_t - 4. * a_t * c - 4. * a * c_t, num = -b + PetscSqrtScalar(disc), /* choose positive sign */
107:     num_t = -b_t + 0.5 / PetscSqrtScalar(disc) * disc_t, den = 2. * a, den_t = 2. * a_t;
108:   *x   = num / den;
109:   *x_t = (num_t * den - num * den_t) / PetscSqr(den);
110: }

112: /* The primary model presented in the paper */
113: static void RDMaterialEnergy_Saha(RD rd, const RDNode *n, PetscScalar *inEm, RDNode *dEm)
114: {
115:   PetscScalar Em, alpha, alpha_t, T = n->T, T_t = 1., chi = rd->I_H / (rd->k * T), chi_t = -chi / T * T_t, a = 1., a_t = 0, b = 4. * rd->m_p / rd->rho * PetscPowScalarReal(2. * PETSC_PI * rd->m_e * rd->I_H / PetscSqr(rd->h), 1.5) * PetscExpScalar(-chi) * PetscPowScalarReal(chi, 1.5), /* Eq 7 */
116:     b_t = -b * chi_t + 1.5 * b / chi * chi_t, c = -b, c_t = -b_t;
117:   QuadraticSolve(a, a_t, b, b_t, c, c_t, &alpha, &alpha_t);      /* Solve Eq 7 for alpha */
118:   Em = rd->k * T / rd->m_p * (1.5 * (1. + alpha) + alpha * chi); /* Eq 6 */
119:   if (inEm) *inEm = Em;
120:   if (dEm) {
121:     dEm->E = 0;
122:     dEm->T = Em / T * T_t + rd->k * T / rd->m_p * (1.5 * alpha_t + alpha_t * chi + alpha * chi_t);
123:   }
124: }
125: /* Reduced ionization model, Eq 30 */
126: static void RDMaterialEnergy_Reduced(RD rd, const RDNode *n, PetscScalar *Em, RDNode *dEm)
127: {
128:   PetscScalar alpha, alpha_t, T = n->T, T_t = 1., chi = -0.3 / T, chi_t = -chi / T * T_t, a = 1., a_t = 0., b = PetscExpScalar(chi), b_t = b * chi_t, c = -b, c_t = -b_t;
129:   QuadraticSolve(a, a_t, b, b_t, c, c_t, &alpha, &alpha_t);
130:   if (Em) *Em = (1. + alpha) * T + 0.3 * alpha;
131:   if (dEm) {
132:     dEm->E = 0;
133:     dEm->T = alpha_t * T + (1. + alpha) * T_t + 0.3 * alpha_t;
134:   }
135: }

137: /* Eq 5 */
138: static void RDSigma_R(RD rd, RDNode *n, PetscScalar *sigma_R, RDNode *dsigma_R)
139: {
140:   *sigma_R    = rd->K_R * rd->rho * PetscPowScalar(n->T, -rd->gamma);
141:   dsigma_R->E = 0;
142:   dsigma_R->T = -rd->gamma * (*sigma_R) / n->T;
143: }

145: /* Eq 4 */
146: static void RDDiffusionCoefficient(RD rd, PetscBool limit, RDNode *n, RDNode *nx, PetscScalar *D_R, RDNode *dD_R, RDNode *dxD_R)
147: {
148:   PetscScalar sigma_R, denom;
149:   RDNode      dsigma_R, ddenom, dxdenom;

151:   RDSigma_R(rd, n, &sigma_R, &dsigma_R);
152:   denom     = 3. * rd->rho * sigma_R + (int)limit * PetscAbsScalar(nx->E) / n->E;
153:   ddenom.E  = -(int)limit * PetscAbsScalar(nx->E) / PetscSqr(n->E);
154:   ddenom.T  = 3. * rd->rho * dsigma_R.T;
155:   dxdenom.E = (int)limit * (PetscRealPart(nx->E) < 0 ? -1. : 1.) / n->E;
156:   dxdenom.T = 0;
157:   *D_R      = rd->c / denom;
158:   if (dD_R) {
159:     dD_R->E = -rd->c / PetscSqr(denom) * ddenom.E;
160:     dD_R->T = -rd->c / PetscSqr(denom) * ddenom.T;
161:   }
162:   if (dxD_R) {
163:     dxD_R->E = -rd->c / PetscSqr(denom) * dxdenom.E;
164:     dxD_R->T = -rd->c / PetscSqr(denom) * dxdenom.T;
165:   }
166: }

168: static PetscErrorCode RDStateView(RD rd, Vec X, Vec Xdot, Vec F)
169: {
170:   DMDALocalInfo info;
171:   PetscInt      i;
172:   const RDNode *x, *xdot, *f;
173:   MPI_Comm      comm;

176:   PetscObjectGetComm((PetscObject)rd->da, &comm);
177:   DMDAGetLocalInfo(rd->da, &info);
178:   DMDAVecGetArrayRead(rd->da, X, (void *)&x);
179:   DMDAVecGetArrayRead(rd->da, Xdot, (void *)&xdot);
180:   DMDAVecGetArrayRead(rd->da, F, (void *)&f);
181:   for (i = info.xs; i < info.xs + info.xm; i++) {
182:     PetscCall(PetscSynchronizedPrintf(comm, "x[%" PetscInt_FMT "] (%10.2G,%10.2G) (%10.2G,%10.2G) (%10.2G,%10.2G)\n", i, (double)PetscRealPart(x[i].E), (double)PetscRealPart(x[i].T), (double)PetscRealPart(xdot[i].E), (double)PetscRealPart(xdot[i].T),
183:                                       (double)PetscRealPart(f[i].E), (double)PetscRealPart(f[i].T)));
184:   }
185:   DMDAVecRestoreArrayRead(rd->da, X, (void *)&x);
186:   DMDAVecRestoreArrayRead(rd->da, Xdot, (void *)&xdot);
187:   DMDAVecRestoreArrayRead(rd->da, F, (void *)&f);
188:   PetscSynchronizedFlush(comm, PETSC_STDOUT);
189:   return 0;
190: }

192: static PetscScalar RDRadiation(RD rd, const RDNode *n, RDNode *dn)
193: {
194:   PetscScalar sigma_p = rd->K_p * rd->rho * PetscPowScalar(n->T, -rd->beta), sigma_p_T = -rd->beta * sigma_p / n->T, tmp = 4. * rd->sigma_b * PetscSqr(PetscSqr(n->T)) / rd->c - n->E, tmp_E = -1., tmp_T = 4. * rd->sigma_b * 4 * n->T * (PetscSqr(n->T)) / rd->c, rad = sigma_p * rd->c * rd->rho * tmp, rad_E = sigma_p * rd->c * rd->rho * tmp_E, rad_T = rd->c * rd->rho * (sigma_p_T * tmp + sigma_p * tmp_T);
195:   if (dn) {
196:     dn->E = rad_E;
197:     dn->T = rad_T;
198:   }
199:   return rad;
200: }

202: static PetscScalar RDDiffusion(RD rd, PetscReal hx, const RDNode x[], PetscInt i, RDNode d[])
203: {
204:   PetscReal   ihx = 1. / hx;
205:   RDNode      n_L, nx_L, n_R, nx_R, dD_L, dxD_L, dD_R, dxD_R, dfluxL[2], dfluxR[2];
206:   PetscScalar D_L, D_R, fluxL, fluxR;

208:   n_L.E  = 0.5 * (x[i - 1].E + x[i].E);
209:   n_L.T  = 0.5 * (x[i - 1].T + x[i].T);
210:   nx_L.E = (x[i].E - x[i - 1].E) / hx;
211:   nx_L.T = (x[i].T - x[i - 1].T) / hx;
212:   RDDiffusionCoefficient(rd, PETSC_TRUE, &n_L, &nx_L, &D_L, &dD_L, &dxD_L);
213:   fluxL       = D_L * nx_L.E;
214:   dfluxL[0].E = -ihx * D_L + (0.5 * dD_L.E - ihx * dxD_L.E) * nx_L.E;
215:   dfluxL[1].E = +ihx * D_L + (0.5 * dD_L.E + ihx * dxD_L.E) * nx_L.E;
216:   dfluxL[0].T = (0.5 * dD_L.T - ihx * dxD_L.T) * nx_L.E;
217:   dfluxL[1].T = (0.5 * dD_L.T + ihx * dxD_L.T) * nx_L.E;

219:   n_R.E  = 0.5 * (x[i].E + x[i + 1].E);
220:   n_R.T  = 0.5 * (x[i].T + x[i + 1].T);
221:   nx_R.E = (x[i + 1].E - x[i].E) / hx;
222:   nx_R.T = (x[i + 1].T - x[i].T) / hx;
223:   RDDiffusionCoefficient(rd, PETSC_TRUE, &n_R, &nx_R, &D_R, &dD_R, &dxD_R);
224:   fluxR       = D_R * nx_R.E;
225:   dfluxR[0].E = -ihx * D_R + (0.5 * dD_R.E - ihx * dxD_R.E) * nx_R.E;
226:   dfluxR[1].E = +ihx * D_R + (0.5 * dD_R.E + ihx * dxD_R.E) * nx_R.E;
227:   dfluxR[0].T = (0.5 * dD_R.T - ihx * dxD_R.T) * nx_R.E;
228:   dfluxR[1].T = (0.5 * dD_R.T + ihx * dxD_R.T) * nx_R.E;

230:   if (d) {
231:     d[0].E = -ihx * dfluxL[0].E;
232:     d[0].T = -ihx * dfluxL[0].T;
233:     d[1].E = ihx * (dfluxR[0].E - dfluxL[1].E);
234:     d[1].T = ihx * (dfluxR[0].T - dfluxL[1].T);
235:     d[2].E = ihx * dfluxR[1].E;
236:     d[2].T = ihx * dfluxR[1].T;
237:   }
238:   return ihx * (fluxR - fluxL);
239: }

241: static PetscErrorCode RDGetLocalArrays(RD rd, TS ts, Vec X, Vec Xdot, PetscReal *Theta, PetscReal *dt, Vec *X0loc, RDNode **x0, Vec *Xloc, RDNode **x, Vec *Xloc_t, RDNode **xdot)
242: {
243:   PetscBool istheta;

246:   DMGetLocalVector(rd->da, X0loc);
247:   DMGetLocalVector(rd->da, Xloc);
248:   DMGetLocalVector(rd->da, Xloc_t);

250:   DMGlobalToLocalBegin(rd->da, X, INSERT_VALUES, *Xloc);
251:   DMGlobalToLocalEnd(rd->da, X, INSERT_VALUES, *Xloc);
252:   DMGlobalToLocalBegin(rd->da, Xdot, INSERT_VALUES, *Xloc_t);
253:   DMGlobalToLocalEnd(rd->da, Xdot, INSERT_VALUES, *Xloc_t);

255:   /*
256:     The following is a hack to subvert TSTHETA which is like an implicit midpoint method to behave more like a trapezoid
257:     rule.  These methods have equivalent linear stability, but the nonlinear stability is somewhat different.  The
258:     radiation system is inconvenient to write in explicit form because the ionization model is "on the left".
259:    */
260:   PetscObjectTypeCompare((PetscObject)ts, TSTHETA, &istheta);
261:   if (istheta && rd->endpoint) {
262:     TSThetaGetTheta(ts, Theta);
263:   } else *Theta = 1.;

265:   TSGetTimeStep(ts, dt);
266:   VecWAXPY(*X0loc, -(*Theta) * (*dt), *Xloc_t, *Xloc); /* back out the value at the start of this step */
267:   if (rd->endpoint) { VecWAXPY(*Xloc, *dt, *Xloc_t, *X0loc); /* move the abscissa to the end of the step */ }

269:   DMDAVecGetArray(rd->da, *X0loc, x0);
270:   DMDAVecGetArray(rd->da, *Xloc, x);
271:   DMDAVecGetArray(rd->da, *Xloc_t, xdot);
272:   return 0;
273: }

275: static PetscErrorCode RDRestoreLocalArrays(RD rd, Vec *X0loc, RDNode **x0, Vec *Xloc, RDNode **x, Vec *Xloc_t, RDNode **xdot)
276: {
278:   DMDAVecRestoreArray(rd->da, *X0loc, x0);
279:   DMDAVecRestoreArray(rd->da, *Xloc, x);
280:   DMDAVecRestoreArray(rd->da, *Xloc_t, xdot);
281:   DMRestoreLocalVector(rd->da, X0loc);
282:   DMRestoreLocalVector(rd->da, Xloc);
283:   DMRestoreLocalVector(rd->da, Xloc_t);
284:   return 0;
285: }

287: static PetscErrorCode PETSC_UNUSED RDCheckDomain_Private(RD rd, TS ts, Vec X, PetscBool *in)
288: {
289:   PetscInt  minloc;
290:   PetscReal min;

293:   VecMin(X, &minloc, &min);
294:   if (min < 0) {
295:     SNES snes;
296:     *in = PETSC_FALSE;
297:     TSGetSNES(ts, &snes);
298:     SNESSetFunctionDomainError(snes);
299:     PetscInfo(ts, "Domain violation at %" PetscInt_FMT " field %" PetscInt_FMT " value %g\n", minloc / 2, minloc % 2, (double)min);
300:   } else *in = PETSC_TRUE;
301:   return 0;
302: }

304: /* Energy and temperature must remain positive */
305: #define RDCheckDomain(rd, ts, X) \
306:   do { \
307:     PetscBool _in; \
308:     RDCheckDomain_Private(rd, ts, X, &_in); \
309:     if (!_in) return 0; \
310:   } while (0)

312: static PetscErrorCode RDIFunction_FD(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
313: {
314:   RD            rd = (RD)ctx;
315:   RDNode       *x, *x0, *xdot, *f;
316:   Vec           X0loc, Xloc, Xloc_t;
317:   PetscReal     hx, Theta, dt;
318:   DMDALocalInfo info;
319:   PetscInt      i;

322:   RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot);
323:   DMDAVecGetArray(rd->da, F, &f);
324:   DMDAGetLocalInfo(rd->da, &info);
325:   VecZeroEntries(F);

327:   hx = rd->L / (info.mx - 1);

329:   for (i = info.xs; i < info.xs + info.xm; i++) {
330:     PetscReal   rho = rd->rho;
331:     PetscScalar Em_t, rad;

333:     rad = (1. - Theta) * RDRadiation(rd, &x0[i], 0) + Theta * RDRadiation(rd, &x[i], 0);
334:     if (rd->endpoint) {
335:       PetscScalar Em0, Em1;
336:       RDMaterialEnergy(rd, &x0[i], &Em0, NULL);
337:       RDMaterialEnergy(rd, &x[i], &Em1, NULL);
338:       Em_t = (Em1 - Em0) / dt;
339:     } else {
340:       RDNode dEm;
341:       RDMaterialEnergy(rd, &x[i], NULL, &dEm);
342:       Em_t = dEm.E * xdot[i].E + dEm.T * xdot[i].T;
343:     }
344:     /* Residuals are multiplied by the volume element (hx).  */
345:     /* The temperature equation does not have boundary conditions */
346:     f[i].T = hx * (rho * Em_t + rad);

348:     if (i == 0) { /* Left boundary condition */
349:       PetscScalar D_R, bcTheta = rd->bcmidpoint ? Theta : 1.;
350:       RDNode      n, nx;

352:       n.E  = (1. - bcTheta) * x0[0].E + bcTheta * x[0].E;
353:       n.T  = (1. - bcTheta) * x0[0].T + bcTheta * x[0].T;
354:       nx.E = ((1. - bcTheta) * (x0[1].E - x0[0].E) + bcTheta * (x[1].E - x[0].E)) / hx;
355:       nx.T = ((1. - bcTheta) * (x0[1].T - x0[0].T) + bcTheta * (x[1].T - x[0].T)) / hx;
356:       switch (rd->leftbc) {
357:       case BC_ROBIN:
358:         RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D_R, 0, 0);
359:         f[0].E = hx * (n.E - 2. * D_R * nx.E - rd->Eapplied);
360:         break;
361:       case BC_NEUMANN:
362:         f[0].E = x[1].E - x[0].E;
363:         break;
364:       default:
365:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT, rd->initial);
366:       }
367:     } else if (i == info.mx - 1) {  /* Right boundary */
368:       f[i].E = x[i].E - x[i - 1].E; /* Homogeneous Neumann */
369:     } else {
370:       PetscScalar diff = (1. - Theta) * RDDiffusion(rd, hx, x0, i, 0) + Theta * RDDiffusion(rd, hx, x, i, 0);
371:       f[i].E           = hx * (xdot[i].E - diff - rad);
372:     }
373:   }
374:   RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot);
375:   DMDAVecRestoreArray(rd->da, F, &f);
376:   if (rd->monitor_residual) RDStateView(rd, X, Xdot, F);
377:   return 0;
378: }

380: static PetscErrorCode RDIJacobian_FD(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ctx)
381: {
382:   RD            rd = (RD)ctx;
383:   RDNode       *x, *x0, *xdot;
384:   Vec           X0loc, Xloc, Xloc_t;
385:   PetscReal     hx, Theta, dt;
386:   DMDALocalInfo info;
387:   PetscInt      i;

390:   RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot);
391:   DMDAGetLocalInfo(rd->da, &info);
392:   hx = rd->L / (info.mx - 1);
393:   MatZeroEntries(B);

395:   for (i = info.xs; i < info.xs + info.xm; i++) {
396:     PetscInt                  col[3];
397:     PetscReal                 rho = rd->rho;
398:     PetscScalar /*Em_t,rad,*/ K[2][6];
399:     RDNode                    dEm_t, drad;

401:     /*rad = (1.-Theta)* */ RDRadiation(rd, &x0[i], 0); /* + Theta* */
402:     RDRadiation(rd, &x[i], &drad);

404:     if (rd->endpoint) {
405:       PetscScalar Em0, Em1;
406:       RDNode      dEm1;
407:       RDMaterialEnergy(rd, &x0[i], &Em0, NULL);
408:       RDMaterialEnergy(rd, &x[i], &Em1, &dEm1);
409:       /*Em_t = (Em1 - Em0) / (Theta*dt);*/
410:       dEm_t.E = dEm1.E / (Theta * dt);
411:       dEm_t.T = dEm1.T / (Theta * dt);
412:     } else {
413:       const PetscScalar epsilon = x[i].T * PETSC_SQRT_MACHINE_EPSILON;
414:       RDNode            n1;
415:       RDNode            dEm, dEm1;
416:       PetscScalar       Em_TT;

418:       n1.E = x[i].E;
419:       n1.T = x[i].T + epsilon;
420:       RDMaterialEnergy(rd, &x[i], NULL, &dEm);
421:       RDMaterialEnergy(rd, &n1, NULL, &dEm1);
422:       /* The Jacobian needs another derivative.  We finite difference here instead of
423:        * propagating second derivatives through the ionization model. */
424:       Em_TT = (dEm1.T - dEm.T) / epsilon;
425:       /*Em_t = dEm.E * xdot[i].E + dEm.T * xdot[i].T;*/
426:       dEm_t.E = dEm.E * a;
427:       dEm_t.T = dEm.T * a + Em_TT * xdot[i].T;
428:     }

430:     PetscMemzero(K, sizeof(K));
431:     /* Residuals are multiplied by the volume element (hx).  */
432:     if (i == 0) {
433:       PetscScalar D, bcTheta = rd->bcmidpoint ? Theta : 1.;
434:       RDNode      n, nx;
435:       RDNode      dD, dxD;

437:       n.E  = (1. - bcTheta) * x0[0].E + bcTheta * x[0].E;
438:       n.T  = (1. - bcTheta) * x0[0].T + bcTheta * x[0].T;
439:       nx.E = ((1. - bcTheta) * (x0[1].E - x0[0].E) + bcTheta * (x[1].E - x[0].E)) / hx;
440:       nx.T = ((1. - bcTheta) * (x0[1].T - x0[0].T) + bcTheta * (x[1].T - x[0].T)) / hx;
441:       switch (rd->leftbc) {
442:       case BC_ROBIN:
443:         RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, &dD, &dxD);
444:         K[0][1 * 2 + 0] = (bcTheta / Theta) * hx * (1. - 2. * D * (-1. / hx) - 2. * nx.E * dD.E + 2. * nx.E * dxD.E / hx);
445:         K[0][1 * 2 + 1] = (bcTheta / Theta) * hx * (-2. * nx.E * dD.T);
446:         K[0][2 * 2 + 0] = (bcTheta / Theta) * hx * (-2. * D * (1. / hx) - 2. * nx.E * dD.E - 2. * nx.E * dxD.E / hx);
447:         break;
448:       case BC_NEUMANN:
449:         K[0][1 * 2 + 0] = -1. / Theta;
450:         K[0][2 * 2 + 0] = 1. / Theta;
451:         break;
452:       default:
453:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT, rd->initial);
454:       }
455:     } else if (i == info.mx - 1) {
456:       K[0][0 * 2 + 0] = -1. / Theta;
457:       K[0][1 * 2 + 0] = 1. / Theta;
458:     } else {
459:       /*PetscScalar diff;*/
460:       RDNode ddiff[3];
461:       /*diff = (1.-Theta)*RDDiffusion(rd,hx,x0,i,0) + Theta* */ RDDiffusion(rd, hx, x, i, ddiff);
462:       K[0][0 * 2 + 0] = -hx * ddiff[0].E;
463:       K[0][0 * 2 + 1] = -hx * ddiff[0].T;
464:       K[0][1 * 2 + 0] = hx * (a - ddiff[1].E - drad.E);
465:       K[0][1 * 2 + 1] = hx * (-ddiff[1].T - drad.T);
466:       K[0][2 * 2 + 0] = -hx * ddiff[2].E;
467:       K[0][2 * 2 + 1] = -hx * ddiff[2].T;
468:     }

470:     K[1][1 * 2 + 0] = hx * (rho * dEm_t.E + drad.E);
471:     K[1][1 * 2 + 1] = hx * (rho * dEm_t.T + drad.T);

473:     col[0] = i - 1;
474:     col[1] = i;
475:     col[2] = i + 1 < info.mx ? i + 1 : -1;
476:     MatSetValuesBlocked(B, 1, &i, 3, col, &K[0][0], INSERT_VALUES);
477:   }
478:   RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot);
479:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
480:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
481:   if (A != B) {
482:     MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
483:     MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
484:   }
485:   return 0;
486: }

488: /* Evaluate interpolants and derivatives at a select quadrature point */
489: static void RDEvaluate(PetscReal interp[][2], PetscReal deriv[][2], PetscInt q, const RDNode x[], PetscInt i, RDNode *n, RDNode *nx)
490: {
491:   PetscInt j;
492:   n->E  = 0;
493:   n->T  = 0;
494:   nx->E = 0;
495:   nx->T = 0;
496:   for (j = 0; j < 2; j++) {
497:     n->E += interp[q][j] * x[i + j].E;
498:     n->T += interp[q][j] * x[i + j].T;
499:     nx->E += deriv[q][j] * x[i + j].E;
500:     nx->T += deriv[q][j] * x[i + j].T;
501:   }
502: }

504: /*
505:  Various quadrature rules.  The nonlinear terms are non-polynomial so no standard quadrature will be exact.
506: */
507: static PetscErrorCode RDGetQuadrature(RD rd, PetscReal hx, PetscInt *nq, PetscReal weight[], PetscReal interp[][2], PetscReal deriv[][2])
508: {
509:   PetscInt         q, j;
510:   const PetscReal *refweight, (*refinterp)[2], (*refderiv)[2];

513:   switch (rd->quadrature) {
514:   case QUADRATURE_GAUSS1: {
515:     static const PetscReal ww[1]    = {1.};
516:     static const PetscReal ii[1][2] = {
517:       {0.5, 0.5}
518:     };
519:     static const PetscReal dd[1][2] = {
520:       {-1., 1.}
521:     };
522:     *nq       = 1;
523:     refweight = ww;
524:     refinterp = ii;
525:     refderiv  = dd;
526:   } break;
527:   case QUADRATURE_GAUSS2: {
528:     static const PetscReal ii[2][2] = {
529:       {0.78867513459481287, 0.21132486540518713},
530:       {0.21132486540518713, 0.78867513459481287}
531:     };
532:     static const PetscReal dd[2][2] = {
533:       {-1., 1.},
534:       {-1., 1.}
535:     };
536:     static const PetscReal ww[2] = {0.5, 0.5};
537:     *nq                          = 2;
538:     refweight                    = ww;
539:     refinterp                    = ii;
540:     refderiv                     = dd;
541:   } break;
542:   case QUADRATURE_GAUSS3: {
543:     static const PetscReal ii[3][2] = {
544:       {0.8872983346207417, 0.1127016653792583},
545:       {0.5,                0.5               },
546:       {0.1127016653792583, 0.8872983346207417}
547:     };
548:     static const PetscReal dd[3][2] = {
549:       {-1, 1},
550:       {-1, 1},
551:       {-1, 1}
552:     };
553:     static const PetscReal ww[3] = {5. / 18, 8. / 18, 5. / 18};
554:     *nq                          = 3;
555:     refweight                    = ww;
556:     refinterp                    = ii;
557:     refderiv                     = dd;
558:   } break;
559:   case QUADRATURE_GAUSS4: {
560:     static const PetscReal ii[][2] = {
561:       {0.93056815579702623,  0.069431844202973658},
562:       {0.66999052179242813,  0.33000947820757187 },
563:       {0.33000947820757187,  0.66999052179242813 },
564:       {0.069431844202973658, 0.93056815579702623 }
565:     };
566:     static const PetscReal dd[][2] = {
567:       {-1, 1},
568:       {-1, 1},
569:       {-1, 1},
570:       {-1, 1}
571:     };
572:     static const PetscReal ww[] = {0.17392742256872692, 0.3260725774312731, 0.3260725774312731, 0.17392742256872692};

574:     *nq       = 4;
575:     refweight = ww;
576:     refinterp = ii;
577:     refderiv  = dd;
578:   } break;
579:   case QUADRATURE_LOBATTO2: {
580:     static const PetscReal ii[2][2] = {
581:       {1., 0.},
582:       {0., 1.}
583:     };
584:     static const PetscReal dd[2][2] = {
585:       {-1., 1.},
586:       {-1., 1.}
587:     };
588:     static const PetscReal ww[2] = {0.5, 0.5};
589:     *nq                          = 2;
590:     refweight                    = ww;
591:     refinterp                    = ii;
592:     refderiv                     = dd;
593:   } break;
594:   case QUADRATURE_LOBATTO3: {
595:     static const PetscReal ii[3][2] = {
596:       {1,   0  },
597:       {0.5, 0.5},
598:       {0,   1  }
599:     };
600:     static const PetscReal dd[3][2] = {
601:       {-1, 1},
602:       {-1, 1},
603:       {-1, 1}
604:     };
605:     static const PetscReal ww[3] = {1. / 6, 4. / 6, 1. / 6};
606:     *nq                          = 3;
607:     refweight                    = ww;
608:     refinterp                    = ii;
609:     refderiv                     = dd;
610:   } break;
611:   default:
612:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown quadrature %d", (int)rd->quadrature);
613:   }

615:   for (q = 0; q < *nq; q++) {
616:     weight[q] = refweight[q] * hx;
617:     for (j = 0; j < 2; j++) {
618:       interp[q][j] = refinterp[q][j];
619:       deriv[q][j]  = refderiv[q][j] / hx;
620:     }
621:   }
622:   return 0;
623: }

625: /*
626:  Finite element version
627: */
628: static PetscErrorCode RDIFunction_FE(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
629: {
630:   RD            rd = (RD)ctx;
631:   RDNode       *x, *x0, *xdot, *f;
632:   Vec           X0loc, Xloc, Xloc_t, Floc;
633:   PetscReal     hx, Theta, dt, weight[5], interp[5][2], deriv[5][2];
634:   DMDALocalInfo info;
635:   PetscInt      i, j, q, nq;

638:   RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot);

640:   DMGetLocalVector(rd->da, &Floc);
641:   VecZeroEntries(Floc);
642:   DMDAVecGetArray(rd->da, Floc, &f);
643:   DMDAGetLocalInfo(rd->da, &info);

645:   /* Set up shape functions and quadrature for elements (assumes a uniform grid) */
646:   hx = rd->L / (info.mx - 1);
647:   RDGetQuadrature(rd, hx, &nq, weight, interp, deriv);

649:   for (i = info.xs; i < PetscMin(info.xs + info.xm, info.mx - 1); i++) {
650:     for (q = 0; q < nq; q++) {
651:       PetscReal   rho = rd->rho;
652:       PetscScalar Em_t, rad, D_R, D0_R;
653:       RDNode      n, n0, nx, n0x, nt, ntx;
654:       RDEvaluate(interp, deriv, q, x, i, &n, &nx);
655:       RDEvaluate(interp, deriv, q, x0, i, &n0, &n0x);
656:       RDEvaluate(interp, deriv, q, xdot, i, &nt, &ntx);

658:       rad = (1. - Theta) * RDRadiation(rd, &n0, 0) + Theta * RDRadiation(rd, &n, 0);
659:       if (rd->endpoint) {
660:         PetscScalar Em0, Em1;
661:         RDMaterialEnergy(rd, &n0, &Em0, NULL);
662:         RDMaterialEnergy(rd, &n, &Em1, NULL);
663:         Em_t = (Em1 - Em0) / dt;
664:       } else {
665:         RDNode dEm;
666:         RDMaterialEnergy(rd, &n, NULL, &dEm);
667:         Em_t = dEm.E * nt.E + dEm.T * nt.T;
668:       }
669:       RDDiffusionCoefficient(rd, PETSC_TRUE, &n0, &n0x, &D0_R, 0, 0);
670:       RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, 0, 0);
671:       for (j = 0; j < 2; j++) {
672:         f[i + j].E += (deriv[q][j] * weight[q] * ((1. - Theta) * D0_R * n0x.E + Theta * D_R * nx.E) + interp[q][j] * weight[q] * (nt.E - rad));
673:         f[i + j].T += interp[q][j] * weight[q] * (rho * Em_t + rad);
674:       }
675:     }
676:   }
677:   if (info.xs == 0) {
678:     switch (rd->leftbc) {
679:     case BC_ROBIN: {
680:       PetscScalar D_R, D_R_bc;
681:       PetscReal   ratio, bcTheta = rd->bcmidpoint ? Theta : 1.;
682:       RDNode      n, nx;

684:       n.E  = (1 - bcTheta) * x0[0].E + bcTheta * x[0].E;
685:       n.T  = (1 - bcTheta) * x0[0].T + bcTheta * x[0].T;
686:       nx.E = (x[1].E - x[0].E) / hx;
687:       nx.T = (x[1].T - x[0].T) / hx;
688:       RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, 0, 0);
689:       RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D_R_bc, 0, 0);
690:       ratio = PetscRealPart(D_R / D_R_bc);
693:       f[0].E += -ratio * 0.5 * (rd->Eapplied - n.E);
694:     } break;
695:     case BC_NEUMANN:
696:       /* homogeneous Neumann is the natural condition */
697:       break;
698:     default:
699:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT, rd->initial);
700:     }
701:   }

703:   RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot);
704:   DMDAVecRestoreArray(rd->da, Floc, &f);
705:   VecZeroEntries(F);
706:   DMLocalToGlobalBegin(rd->da, Floc, ADD_VALUES, F);
707:   DMLocalToGlobalEnd(rd->da, Floc, ADD_VALUES, F);
708:   DMRestoreLocalVector(rd->da, &Floc);

710:   if (rd->monitor_residual) RDStateView(rd, X, Xdot, F);
711:   return 0;
712: }

714: static PetscErrorCode RDIJacobian_FE(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, void *ctx)
715: {
716:   RD            rd = (RD)ctx;
717:   RDNode       *x, *x0, *xdot;
718:   Vec           X0loc, Xloc, Xloc_t;
719:   PetscReal     hx, Theta, dt, weight[5], interp[5][2], deriv[5][2];
720:   DMDALocalInfo info;
721:   PetscInt      i, j, k, q, nq;
722:   PetscScalar   K[4][4];

725:   RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot);
726:   DMDAGetLocalInfo(rd->da, &info);
727:   hx = rd->L / (info.mx - 1);
728:   RDGetQuadrature(rd, hx, &nq, weight, interp, deriv);
729:   MatZeroEntries(B);
730:   for (i = info.xs; i < PetscMin(info.xs + info.xm, info.mx - 1); i++) {
731:     PetscInt rc[2];

733:     rc[0] = i;
734:     rc[1] = i + 1;
735:     PetscMemzero(K, sizeof(K));
736:     for (q = 0; q < nq; q++) {
737:       PetscScalar              D_R;
738:       PETSC_UNUSED PetscScalar rad;
739:       RDNode                   n, nx, nt, ntx, drad, dD_R, dxD_R, dEm;
740:       RDEvaluate(interp, deriv, q, x, i, &n, &nx);
741:       RDEvaluate(interp, deriv, q, xdot, i, &nt, &ntx);
742:       rad = RDRadiation(rd, &n, &drad);
743:       RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, &dD_R, &dxD_R);
744:       RDMaterialEnergy(rd, &n, NULL, &dEm);
745:       for (j = 0; j < 2; j++) {
746:         for (k = 0; k < 2; k++) {
747:           K[j * 2 + 0][k * 2 + 0] += (+interp[q][j] * weight[q] * (a - drad.E) * interp[q][k] + deriv[q][j] * weight[q] * ((D_R + dxD_R.E * nx.E) * deriv[q][k] + dD_R.E * nx.E * interp[q][k]));
748:           K[j * 2 + 0][k * 2 + 1] += (+interp[q][j] * weight[q] * (-drad.T * interp[q][k]) + deriv[q][j] * weight[q] * (dxD_R.T * deriv[q][k] + dD_R.T * interp[q][k]) * nx.E);
749:           K[j * 2 + 1][k * 2 + 0] += interp[q][j] * weight[q] * drad.E * interp[q][k];
750:           K[j * 2 + 1][k * 2 + 1] += interp[q][j] * weight[q] * (a * rd->rho * dEm.T + drad.T) * interp[q][k];
751:         }
752:       }
753:     }
754:     MatSetValuesBlocked(B, 2, rc, 2, rc, &K[0][0], ADD_VALUES);
755:   }
756:   if (info.xs == 0) {
757:     switch (rd->leftbc) {
758:     case BC_ROBIN: {
759:       PetscScalar D_R, D_R_bc;
760:       PetscReal   ratio;
761:       RDNode      n, nx;

763:       n.E  = (1 - Theta) * x0[0].E + Theta * x[0].E;
764:       n.T  = (1 - Theta) * x0[0].T + Theta * x[0].T;
765:       nx.E = (x[1].E - x[0].E) / hx;
766:       nx.T = (x[1].T - x[0].T) / hx;
767:       RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, 0, 0);
768:       RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D_R_bc, 0, 0);
769:       ratio = PetscRealPart(D_R / D_R_bc);
770:       MatSetValue(B, 0, 0, ratio * 0.5, ADD_VALUES);
771:     } break;
772:     case BC_NEUMANN:
773:       /* homogeneous Neumann is the natural condition */
774:       break;
775:     default:
776:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT, rd->initial);
777:     }
778:   }

780:   RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot);
781:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
782:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
783:   if (A != B) {
784:     MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
785:     MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
786:   }
787:   return 0;
788: }

790: /* Temperature that is in equilibrium with the radiation density */
791: static PetscScalar RDRadiationTemperature(RD rd, PetscScalar E)
792: {
793:   return PetscPowScalar(E * rd->c / (4. * rd->sigma_b), 0.25);
794: }

796: static PetscErrorCode RDInitialState(RD rd, Vec X)
797: {
798:   DMDALocalInfo info;
799:   PetscInt      i;
800:   RDNode       *x;

803:   DMDAGetLocalInfo(rd->da, &info);
804:   DMDAVecGetArray(rd->da, X, &x);
805:   for (i = info.xs; i < info.xs + info.xm; i++) {
806:     PetscReal coord = i * rd->L / (info.mx - 1);
807:     switch (rd->initial) {
808:     case 1:
809:       x[i].E = 0.001;
810:       x[i].T = RDRadiationTemperature(rd, x[i].E);
811:       break;
812:     case 2:
813:       x[i].E = 0.001 + 100. * PetscExpReal(-PetscSqr(coord / 0.1));
814:       x[i].T = RDRadiationTemperature(rd, x[i].E);
815:       break;
816:     case 3:
817:       x[i].E = 7.56e-2 * rd->unit.Joule / PetscPowScalarInt(rd->unit.meter, 3);
818:       x[i].T = RDRadiationTemperature(rd, x[i].E);
819:       break;
820:     default:
821:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No initial state %" PetscInt_FMT, rd->initial);
822:     }
823:   }
824:   DMDAVecRestoreArray(rd->da, X, &x);
825:   return 0;
826: }

828: static PetscErrorCode RDView(RD rd, Vec X, PetscViewer viewer)
829: {
830:   Vec             Y;
831:   const RDNode   *x;
832:   PetscScalar    *y;
833:   PetscInt        i, m, M;
834:   const PetscInt *lx;
835:   DM              da;
836:   MPI_Comm        comm;

839:   /*
840:     Create a DMDA (one dof per node, zero stencil width, same layout) to hold Trad
841:     (radiation temperature).  It is not necessary to create a DMDA for this, but this way
842:     output and visualization will have meaningful variable names and correct scales.
843:   */
844:   DMDAGetInfo(rd->da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
845:   DMDAGetOwnershipRanges(rd->da, &lx, 0, 0);
846:   PetscObjectGetComm((PetscObject)rd->da, &comm);
847:   DMDACreate1d(comm, DM_BOUNDARY_NONE, M, 1, 0, lx, &da);
848:   DMSetFromOptions(da);
849:   DMSetUp(da);
850:   DMDASetUniformCoordinates(da, 0., rd->L, 0., 0., 0., 0.);
851:   DMDASetFieldName(da, 0, "T_rad");
852:   DMCreateGlobalVector(da, &Y);

854:   /* Compute the radiation temperature from the solution at each node */
855:   VecGetLocalSize(Y, &m);
856:   VecGetArrayRead(X, (const PetscScalar **)&x);
857:   VecGetArray(Y, &y);
858:   for (i = 0; i < m; i++) y[i] = RDRadiationTemperature(rd, x[i].E);
859:   VecRestoreArrayRead(X, (const PetscScalar **)&x);
860:   VecRestoreArray(Y, &y);

862:   VecView(Y, viewer);
863:   VecDestroy(&Y);
864:   DMDestroy(&da);
865:   return 0;
866: }

868: static PetscErrorCode RDTestDifferentiation(RD rd)
869: {
870:   MPI_Comm    comm;
871:   RDNode      n, nx;
872:   PetscScalar epsilon;

875:   PetscObjectGetComm((PetscObject)rd->da, &comm);
876:   epsilon = 1e-8;
877:   {
878:     RDNode      dEm, fdEm;
879:     PetscScalar T0 = 1000., T1 = T0 * (1. + epsilon), Em0, Em1;
880:     n.E = 1.;
881:     n.T = T0;
882:     rd->MaterialEnergy(rd, &n, &Em0, &dEm);
883:     n.E = 1. + epsilon;
884:     n.T = T0;
885:     rd->MaterialEnergy(rd, &n, &Em1, 0);
886:     fdEm.E = (Em1 - Em0) / epsilon;
887:     n.E    = 1.;
888:     n.T    = T1;
889:     rd->MaterialEnergy(rd, &n, &Em1, 0);
890:     fdEm.T = (Em1 - Em0) / (T0 * epsilon);
891:     PetscCall(PetscPrintf(comm, "dEm {%g,%g}, fdEm {%g,%g}, diff {%g,%g}\n", (double)PetscRealPart(dEm.E), (double)PetscRealPart(dEm.T), (double)PetscRealPart(fdEm.E), (double)PetscRealPart(fdEm.T), (double)PetscRealPart(dEm.E - fdEm.E),
892:                           (double)PetscRealPart(dEm.T - fdEm.T)));
893:   }
894:   {
895:     PetscScalar D0, D;
896:     RDNode      dD, dxD, fdD, fdxD;
897:     n.E  = 1.;
898:     n.T  = 1.;
899:     nx.E = 1.;
900:     n.T  = 1.;
901:     RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D0, &dD, &dxD);
902:     n.E  = 1. + epsilon;
903:     n.T  = 1.;
904:     nx.E = 1.;
905:     n.T  = 1.;
906:     RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0);
907:     fdD.E = (D - D0) / epsilon;
908:     n.E   = 1;
909:     n.T   = 1. + epsilon;
910:     nx.E  = 1.;
911:     n.T   = 1.;
912:     RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0);
913:     fdD.T = (D - D0) / epsilon;
914:     n.E   = 1;
915:     n.T   = 1.;
916:     nx.E  = 1. + epsilon;
917:     n.T   = 1.;
918:     RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0);
919:     fdxD.E = (D - D0) / epsilon;
920:     n.E    = 1;
921:     n.T    = 1.;
922:     nx.E   = 1.;
923:     n.T    = 1. + epsilon;
924:     RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0);
925:     fdxD.T = (D - D0) / epsilon;
926:     PetscCall(PetscPrintf(comm, "dD {%g,%g}, fdD {%g,%g}, diff {%g,%g}\n", (double)PetscRealPart(dD.E), (double)PetscRealPart(dD.T), (double)PetscRealPart(fdD.E), (double)PetscRealPart(fdD.T), (double)PetscRealPart(dD.E - fdD.E),
927:                           (double)PetscRealPart(dD.T - fdD.T)));
928:     PetscCall(PetscPrintf(comm, "dxD {%g,%g}, fdxD {%g,%g}, diffx {%g,%g}\n", (double)PetscRealPart(dxD.E), (double)PetscRealPart(dxD.T), (double)PetscRealPart(fdxD.E), (double)PetscRealPart(fdxD.T), (double)PetscRealPart(dxD.E - fdxD.E),
929:                           (double)PetscRealPart(dxD.T - fdxD.T)));
930:   }
931:   {
932:     PetscInt    i;
933:     PetscReal   hx = 1.;
934:     PetscScalar a0;
935:     RDNode      n0[3], n1[3], d[3], fd[3];

937:     n0[0].E = 1.;
938:     n0[0].T = 1.;
939:     n0[1].E = 5.;
940:     n0[1].T = 3.;
941:     n0[2].E = 4.;
942:     n0[2].T = 2.;
943:     a0      = RDDiffusion(rd, hx, n0, 1, d);
944:     for (i = 0; i < 3; i++) {
945:       PetscMemcpy(n1, n0, sizeof(n0));
946:       n1[i].E += epsilon;
947:       fd[i].E = (RDDiffusion(rd, hx, n1, 1, 0) - a0) / epsilon;
948:       PetscMemcpy(n1, n0, sizeof(n0));
949:       n1[i].T += epsilon;
950:       fd[i].T = (RDDiffusion(rd, hx, n1, 1, 0) - a0) / epsilon;
951:       PetscCall(PetscPrintf(comm, "ddiff[%" PetscInt_FMT "] {%g,%g}, fd {%g %g}, diff {%g,%g}\n", i, (double)PetscRealPart(d[i].E), (double)PetscRealPart(d[i].T), (double)PetscRealPart(fd[i].E), (double)PetscRealPart(fd[i].T),
952:                             (double)PetscRealPart(d[i].E - fd[i].E), (double)PetscRealPart(d[i].T - fd[i].T)));
953:     }
954:   }
955:   {
956:     PetscScalar rad0, rad;
957:     RDNode      drad, fdrad;
958:     n.E     = 1.;
959:     n.T     = 1.;
960:     rad0    = RDRadiation(rd, &n, &drad);
961:     n.E     = 1. + epsilon;
962:     n.T     = 1.;
963:     rad     = RDRadiation(rd, &n, 0);
964:     fdrad.E = (rad - rad0) / epsilon;
965:     n.E     = 1.;
966:     n.T     = 1. + epsilon;
967:     rad     = RDRadiation(rd, &n, 0);
968:     fdrad.T = (rad - rad0) / epsilon;
969:     PetscCall(PetscPrintf(comm, "drad {%g,%g}, fdrad {%g,%g}, diff {%g,%g}\n", (double)PetscRealPart(drad.E), (double)PetscRealPart(drad.T), (double)PetscRealPart(fdrad.E), (double)PetscRealPart(fdrad.T), (double)PetscRealPart(drad.E - drad.E),
970:                           (double)PetscRealPart(drad.T - fdrad.T)));
971:   }
972:   return 0;
973: }

975: static PetscErrorCode RDCreate(MPI_Comm comm, RD *inrd)
976: {
977:   RD        rd;
978:   PetscReal meter = 0, kilogram = 0, second = 0, Kelvin = 0, Joule = 0, Watt = 0;

981:   *inrd = 0;
982:   PetscNew(&rd);

984:   PetscOptionsBegin(comm, NULL, "Options for nonequilibrium radiation-diffusion with RD ionization", NULL);
985:   {
986:     rd->initial = 1;
987:     PetscOptionsInt("-rd_initial", "Initial condition (1=Marshak, 2=Blast, 3=Marshak+)", "", rd->initial, &rd->initial, 0);
988:     switch (rd->initial) {
989:     case 1:
990:     case 2:
991:       rd->unit.kilogram = 1.;
992:       rd->unit.meter    = 1.;
993:       rd->unit.second   = 1.;
994:       rd->unit.Kelvin   = 1.;
995:       break;
996:     case 3:
997:       rd->unit.kilogram = 1.e12;
998:       rd->unit.meter    = 1.;
999:       rd->unit.second   = 1.e9;
1000:       rd->unit.Kelvin   = 1.;
1001:       break;
1002:     default:
1003:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown initial condition %" PetscInt_FMT, rd->initial);
1004:     }
1005:     /* Fundamental units */
1006:     PetscOptionsReal("-rd_unit_meter", "Length of 1 meter in nondimensional units", "", rd->unit.meter, &rd->unit.meter, 0);
1007:     PetscOptionsReal("-rd_unit_kilogram", "Mass of 1 kilogram in nondimensional units", "", rd->unit.kilogram, &rd->unit.kilogram, 0);
1008:     PetscOptionsReal("-rd_unit_second", "Time of a second in nondimensional units", "", rd->unit.second, &rd->unit.second, 0);
1009:     PetscOptionsReal("-rd_unit_Kelvin", "Temperature of a Kelvin in nondimensional units", "", rd->unit.Kelvin, &rd->unit.Kelvin, 0);
1010:     /* Derived units */
1011:     rd->unit.Joule = rd->unit.kilogram * PetscSqr(rd->unit.meter / rd->unit.second);
1012:     rd->unit.Watt  = rd->unit.Joule / rd->unit.second;
1013:     /* Local aliases */
1014:     meter    = rd->unit.meter;
1015:     kilogram = rd->unit.kilogram;
1016:     second   = rd->unit.second;
1017:     Kelvin   = rd->unit.Kelvin;
1018:     Joule    = rd->unit.Joule;
1019:     Watt     = rd->unit.Watt;

1021:     PetscOptionsBool("-rd_monitor_residual", "Display residuals every time they are evaluated", "", rd->monitor_residual, &rd->monitor_residual, NULL);
1022:     PetscOptionsEnum("-rd_discretization", "Discretization type", "", DiscretizationTypes, (PetscEnum)rd->discretization, (PetscEnum *)&rd->discretization, NULL);
1023:     if (rd->discretization == DISCRETIZATION_FE) {
1024:       rd->quadrature = QUADRATURE_GAUSS2;
1025:       PetscOptionsEnum("-rd_quadrature", "Finite element quadrature", "", QuadratureTypes, (PetscEnum)rd->quadrature, (PetscEnum *)&rd->quadrature, NULL);
1026:     }
1027:     PetscOptionsEnum("-rd_jacobian", "Type of finite difference Jacobian", "", JacobianTypes, (PetscEnum)rd->jacobian, (PetscEnum *)&rd->jacobian, NULL);
1028:     switch (rd->initial) {
1029:     case 1:
1030:       rd->leftbc     = BC_ROBIN;
1031:       rd->Eapplied   = 4 * rd->unit.Joule / PetscPowRealInt(rd->unit.meter, 3);
1032:       rd->L          = 1. * rd->unit.meter;
1033:       rd->beta       = 3.0;
1034:       rd->gamma      = 3.0;
1035:       rd->final_time = 3 * second;
1036:       break;
1037:     case 2:
1038:       rd->leftbc     = BC_NEUMANN;
1039:       rd->Eapplied   = 0.;
1040:       rd->L          = 1. * rd->unit.meter;
1041:       rd->beta       = 3.0;
1042:       rd->gamma      = 3.0;
1043:       rd->final_time = 1 * second;
1044:       break;
1045:     case 3:
1046:       rd->leftbc     = BC_ROBIN;
1047:       rd->Eapplied   = 7.503e6 * rd->unit.Joule / PetscPowRealInt(rd->unit.meter, 3);
1048:       rd->L          = 5. * rd->unit.meter;
1049:       rd->beta       = 3.5;
1050:       rd->gamma      = 3.5;
1051:       rd->final_time = 20e-9 * second;
1052:       break;
1053:     default:
1054:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Initial %" PetscInt_FMT, rd->initial);
1055:     }
1056:     PetscOptionsEnum("-rd_leftbc", "Left boundary condition", "", BCTypes, (PetscEnum)rd->leftbc, (PetscEnum *)&rd->leftbc, NULL);
1057:     PetscOptionsReal("-rd_E_applied", "Radiation flux at left end of domain", "", rd->Eapplied, &rd->Eapplied, NULL);
1058:     PetscOptionsReal("-rd_beta", "Thermal exponent for photon absorption", "", rd->beta, &rd->beta, NULL);
1059:     PetscOptionsReal("-rd_gamma", "Thermal exponent for diffusion coefficient", "", rd->gamma, &rd->gamma, NULL);
1060:     PetscOptionsBool("-rd_view_draw", "Draw final solution", "", rd->view_draw, &rd->view_draw, NULL);
1061:     PetscOptionsBool("-rd_endpoint", "Discretize using endpoints (like trapezoid rule) instead of midpoint", "", rd->endpoint, &rd->endpoint, NULL);
1062:     PetscOptionsBool("-rd_bcmidpoint", "Impose the boundary condition at the midpoint (Theta) of the interval", "", rd->bcmidpoint, &rd->bcmidpoint, NULL);
1063:     PetscOptionsBool("-rd_bclimit", "Limit diffusion coefficient in definition of Robin boundary condition", "", rd->bclimit, &rd->bclimit, NULL);
1064:     PetscOptionsBool("-rd_test_diff", "Test differentiation in constitutive relations", "", rd->test_diff, &rd->test_diff, NULL);
1065:     PetscOptionsString("-rd_view_binary", "File name to hold final solution", "", rd->view_binary, rd->view_binary, sizeof(rd->view_binary), NULL);
1066:   }
1067:   PetscOptionsEnd();

1069:   switch (rd->initial) {
1070:   case 1:
1071:   case 2:
1072:     rd->rho            = 1.;
1073:     rd->c              = 1.;
1074:     rd->K_R            = 1.;
1075:     rd->K_p            = 1.;
1076:     rd->sigma_b        = 0.25;
1077:     rd->MaterialEnergy = RDMaterialEnergy_Reduced;
1078:     break;
1079:   case 3:
1080:     /* Table 2 */
1081:     rd->rho            = 1.17e-3 * kilogram / (meter * meter * meter);                                                    /* density */
1082:     rd->K_R            = 7.44e18 * PetscPowRealInt(meter, 5) * PetscPowReal(Kelvin, 3.5) * PetscPowRealInt(kilogram, -2); /*  */
1083:     rd->K_p            = 2.33e20 * PetscPowRealInt(meter, 5) * PetscPowReal(Kelvin, 3.5) * PetscPowRealInt(kilogram, -2); /*  */
1084:     rd->I_H            = 2.179e-18 * Joule;                                                                               /* Hydrogen ionization potential */
1085:     rd->m_p            = 1.673e-27 * kilogram;                                                                            /* proton mass */
1086:     rd->m_e            = 9.109e-31 * kilogram;                                                                            /* electron mass */
1087:     rd->h              = 6.626e-34 * Joule * second;                                                                      /* Planck's constant */
1088:     rd->k              = 1.381e-23 * Joule / Kelvin;                                                                      /* Boltzman constant */
1089:     rd->c              = 3.00e8 * meter / second;                                                                         /* speed of light */
1090:     rd->sigma_b        = 5.67e-8 * Watt * PetscPowRealInt(meter, -2) * PetscPowRealInt(Kelvin, -4);                       /* Stefan-Boltzman constant */
1091:     rd->MaterialEnergy = RDMaterialEnergy_Saha;
1092:     break;
1093:   }

1095:   DMDACreate1d(comm, DM_BOUNDARY_NONE, 20, sizeof(RDNode) / sizeof(PetscScalar), 1, NULL, &rd->da);
1096:   DMSetFromOptions(rd->da);
1097:   DMSetUp(rd->da);
1098:   DMDASetFieldName(rd->da, 0, "E");
1099:   DMDASetFieldName(rd->da, 1, "T");
1100:   DMDASetUniformCoordinates(rd->da, 0., 1., 0., 0., 0., 0.);

1102:   *inrd = rd;
1103:   return 0;
1104: }

1106: int main(int argc, char *argv[])
1107: {
1108:   RD        rd;
1109:   TS        ts;
1110:   SNES      snes;
1111:   Vec       X;
1112:   Mat       A, B;
1113:   PetscInt  steps;
1114:   PetscReal ftime;

1117:   PetscInitialize(&argc, &argv, 0, NULL);
1118:   RDCreate(PETSC_COMM_WORLD, &rd);
1119:   DMCreateGlobalVector(rd->da, &X);
1120:   DMSetMatType(rd->da, MATAIJ);
1121:   DMCreateMatrix(rd->da, &B);
1122:   RDInitialState(rd, X);

1124:   TSCreate(PETSC_COMM_WORLD, &ts);
1125:   TSSetProblemType(ts, TS_NONLINEAR);
1126:   TSSetType(ts, TSTHETA);
1127:   TSSetDM(ts, rd->da);
1128:   switch (rd->discretization) {
1129:   case DISCRETIZATION_FD:
1130:     TSSetIFunction(ts, NULL, RDIFunction_FD, rd);
1131:     if (rd->jacobian == JACOBIAN_ANALYTIC) TSSetIJacobian(ts, B, B, RDIJacobian_FD, rd);
1132:     break;
1133:   case DISCRETIZATION_FE:
1134:     TSSetIFunction(ts, NULL, RDIFunction_FE, rd);
1135:     if (rd->jacobian == JACOBIAN_ANALYTIC) TSSetIJacobian(ts, B, B, RDIJacobian_FE, rd);
1136:     break;
1137:   }
1138:   TSSetMaxTime(ts, rd->final_time);
1139:   TSSetTimeStep(ts, 1e-3);
1140:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1141:   TSSetFromOptions(ts);

1143:   A = B;
1144:   TSGetSNES(ts, &snes);
1145:   switch (rd->jacobian) {
1146:   case JACOBIAN_ANALYTIC:
1147:     break;
1148:   case JACOBIAN_MATRIXFREE:
1149:     break;
1150:   case JACOBIAN_FD_COLORING: {
1151:     SNESSetJacobian(snes, A, B, SNESComputeJacobianDefaultColor, 0);
1152:   } break;
1153:   case JACOBIAN_FD_FULL:
1154:     SNESSetJacobian(snes, A, B, SNESComputeJacobianDefault, ts);
1155:     break;
1156:   }

1158:   if (rd->test_diff) RDTestDifferentiation(rd);
1159:   TSSolve(ts, X);
1160:   TSGetSolveTime(ts, &ftime);
1161:   TSGetStepNumber(ts, &steps);
1162:   PetscPrintf(PETSC_COMM_WORLD, "Steps %" PetscInt_FMT "  final time %g\n", steps, (double)ftime);
1163:   if (rd->view_draw) RDView(rd, X, PETSC_VIEWER_DRAW_WORLD);
1164:   if (rd->view_binary[0]) {
1165:     PetscViewer viewer;
1166:     PetscViewerBinaryOpen(PETSC_COMM_WORLD, rd->view_binary, FILE_MODE_WRITE, &viewer);
1167:     RDView(rd, X, viewer);
1168:     PetscViewerDestroy(&viewer);
1169:   }
1170:   VecDestroy(&X);
1171:   MatDestroy(&B);
1172:   RDDestroy(&rd);
1173:   TSDestroy(&ts);
1174:   PetscFinalize();
1175:   return 0;
1176: }
1177: /*TEST

1179:     test:
1180:       args: -da_grid_x 20 -rd_initial 1 -rd_discretization fd -rd_jacobian fd_coloring -rd_endpoint -ts_max_time 1 -ts_dt 2e-1 -ts_theta_initial_guess_extrapolate 0 -ts_monitor -snes_monitor_short -ksp_monitor_short
1181:       requires: !single

1183:     test:
1184:       suffix: 2
1185:       args: -da_grid_x 20 -rd_initial 1 -rd_discretization fe -rd_quadrature lobatto2 -rd_jacobian fd_coloring -rd_endpoint -ts_max_time 1 -ts_dt 2e-1 -ts_theta_initial_guess_extrapolate 0 -ts_monitor -snes_monitor_short -ksp_monitor_short
1186:       requires: !single

1188:     test:
1189:       suffix: 3
1190:       nsize: 2
1191:       args: -da_grid_x 20 -rd_initial 1 -rd_discretization fd -rd_jacobian analytic -rd_endpoint -ts_max_time 3 -ts_dt 1e-1 -ts_theta_initial_guess_extrapolate 0 -ts_monitor -snes_monitor_short -ksp_monitor_short
1192:       requires: !single

1194: TEST*/