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: {
84: PetscFunctionBeginUser;
85: PetscCall(DMDestroy(&(*rd)->da));
86: PetscCall(PetscFree(*rd));
87: PetscFunctionReturn(PETSC_SUCCESS);
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;
175: PetscFunctionBeginUser;
176: PetscCall(PetscObjectGetComm((PetscObject)rd->da, &comm));
177: PetscCall(DMDAGetLocalInfo(rd->da, &info));
178: PetscCall(DMDAVecGetArrayRead(rd->da, X, (void *)&x));
179: PetscCall(DMDAVecGetArrayRead(rd->da, Xdot, (void *)&xdot));
180: PetscCall(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: PetscCall(DMDAVecRestoreArrayRead(rd->da, X, (void *)&x));
186: PetscCall(DMDAVecRestoreArrayRead(rd->da, Xdot, (void *)&xdot));
187: PetscCall(DMDAVecRestoreArrayRead(rd->da, F, (void *)&f));
188: PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
189: PetscFunctionReturn(PETSC_SUCCESS);
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;
245: PetscFunctionBeginUser;
246: PetscCall(DMGetLocalVector(rd->da, X0loc));
247: PetscCall(DMGetLocalVector(rd->da, Xloc));
248: PetscCall(DMGetLocalVector(rd->da, Xloc_t));
250: PetscCall(DMGlobalToLocalBegin(rd->da, X, INSERT_VALUES, *Xloc));
251: PetscCall(DMGlobalToLocalEnd(rd->da, X, INSERT_VALUES, *Xloc));
252: PetscCall(DMGlobalToLocalBegin(rd->da, Xdot, INSERT_VALUES, *Xloc_t));
253: PetscCall(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: PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSTHETA, &istheta));
261: if (istheta && rd->endpoint) {
262: PetscCall(TSThetaGetTheta(ts, Theta));
263: } else *Theta = 1.;
265: PetscCall(TSGetTimeStep(ts, dt));
266: PetscCall(VecWAXPY(*X0loc, -(*Theta) * (*dt), *Xloc_t, *Xloc)); /* back out the value at the start of this step */
267: if (rd->endpoint) { PetscCall(VecWAXPY(*Xloc, *dt, *Xloc_t, *X0loc)); /* move the abscissa to the end of the step */ }
269: PetscCall(DMDAVecGetArray(rd->da, *X0loc, x0));
270: PetscCall(DMDAVecGetArray(rd->da, *Xloc, x));
271: PetscCall(DMDAVecGetArray(rd->da, *Xloc_t, xdot));
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: static PetscErrorCode RDRestoreLocalArrays(RD rd, Vec *X0loc, RDNode **x0, Vec *Xloc, RDNode **x, Vec *Xloc_t, RDNode **xdot)
276: {
277: PetscFunctionBeginUser;
278: PetscCall(DMDAVecRestoreArray(rd->da, *X0loc, x0));
279: PetscCall(DMDAVecRestoreArray(rd->da, *Xloc, x));
280: PetscCall(DMDAVecRestoreArray(rd->da, *Xloc_t, xdot));
281: PetscCall(DMRestoreLocalVector(rd->da, X0loc));
282: PetscCall(DMRestoreLocalVector(rd->da, Xloc));
283: PetscCall(DMRestoreLocalVector(rd->da, Xloc_t));
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: static PetscErrorCode PETSC_UNUSED RDCheckDomain_Private(RD rd, TS ts, Vec X, PetscBool *in)
288: {
289: PetscInt minloc;
290: PetscReal min;
292: PetscFunctionBeginUser;
293: PetscCall(VecMin(X, &minloc, &min));
294: if (min < 0) {
295: SNES snes;
296: *in = PETSC_FALSE;
297: PetscCall(TSGetSNES(ts, &snes));
298: PetscCall(SNESSetFunctionDomainError(snes));
299: PetscCall(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: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: /* Energy and temperature must remain positive */
305: #define RDCheckDomain(rd, ts, X) \
306: do { \
307: PetscBool _in; \
308: PetscCall(RDCheckDomain_Private(rd, ts, X, &_in)); \
309: if (!_in) PetscFunctionReturn(PETSC_SUCCESS); \
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;
321: PetscFunctionBeginUser;
322: PetscCall(RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot));
323: PetscCall(DMDAVecGetArray(rd->da, F, &f));
324: PetscCall(DMDAGetLocalInfo(rd->da, &info));
325: PetscCall(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: PetscCall(RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot));
375: PetscCall(DMDAVecRestoreArray(rd->da, F, &f));
376: if (rd->monitor_residual) PetscCall(RDStateView(rd, X, Xdot, F));
377: PetscFunctionReturn(PETSC_SUCCESS);
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;
389: PetscFunctionBeginUser;
390: PetscCall(RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot));
391: PetscCall(DMDAGetLocalInfo(rd->da, &info));
392: hx = rd->L / (info.mx - 1);
393: PetscCall(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: PetscCall(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: PetscCall(MatSetValuesBlocked(B, 1, &i, 3, col, &K[0][0], INSERT_VALUES));
477: }
478: PetscCall(RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot));
479: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
480: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
481: if (A != B) {
482: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
483: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
484: }
485: PetscFunctionReturn(PETSC_SUCCESS);
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];
512: PetscFunctionBeginUser;
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: PetscFunctionReturn(PETSC_SUCCESS);
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;
637: PetscFunctionBeginUser;
638: PetscCall(RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot));
640: PetscCall(DMGetLocalVector(rd->da, &Floc));
641: PetscCall(VecZeroEntries(Floc));
642: PetscCall(DMDAVecGetArray(rd->da, Floc, &f));
643: PetscCall(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: PetscCall(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);
691: PetscCheck(ratio <= 1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Limited diffusivity is greater than unlimited");
692: PetscCheck(ratio >= 1e-3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Heavily limited diffusivity");
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: PetscCall(RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot));
704: PetscCall(DMDAVecRestoreArray(rd->da, Floc, &f));
705: PetscCall(VecZeroEntries(F));
706: PetscCall(DMLocalToGlobalBegin(rd->da, Floc, ADD_VALUES, F));
707: PetscCall(DMLocalToGlobalEnd(rd->da, Floc, ADD_VALUES, F));
708: PetscCall(DMRestoreLocalVector(rd->da, &Floc));
710: if (rd->monitor_residual) PetscCall(RDStateView(rd, X, Xdot, F));
711: PetscFunctionReturn(PETSC_SUCCESS);
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];
724: PetscFunctionBeginUser;
725: PetscCall(RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot));
726: PetscCall(DMDAGetLocalInfo(rd->da, &info));
727: hx = rd->L / (info.mx - 1);
728: PetscCall(RDGetQuadrature(rd, hx, &nq, weight, interp, deriv));
729: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot));
781: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
782: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
783: if (A != B) {
784: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
785: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
786: }
787: PetscFunctionReturn(PETSC_SUCCESS);
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;
802: PetscFunctionBeginUser;
803: PetscCall(DMDAGetLocalInfo(rd->da, &info));
804: PetscCall(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: PetscCall(DMDAVecRestoreArray(rd->da, X, &x));
825: PetscFunctionReturn(PETSC_SUCCESS);
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;
838: PetscFunctionBeginUser;
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: PetscCall(DMDAGetInfo(rd->da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
845: PetscCall(DMDAGetOwnershipRanges(rd->da, &lx, 0, 0));
846: PetscCall(PetscObjectGetComm((PetscObject)rd->da, &comm));
847: PetscCall(DMDACreate1d(comm, DM_BOUNDARY_NONE, M, 1, 0, lx, &da));
848: PetscCall(DMSetFromOptions(da));
849: PetscCall(DMSetUp(da));
850: PetscCall(DMDASetUniformCoordinates(da, 0., rd->L, 0., 0., 0., 0.));
851: PetscCall(DMDASetFieldName(da, 0, "T_rad"));
852: PetscCall(DMCreateGlobalVector(da, &Y));
854: /* Compute the radiation temperature from the solution at each node */
855: PetscCall(VecGetLocalSize(Y, &m));
856: PetscCall(VecGetArrayRead(X, (const PetscScalar **)&x));
857: PetscCall(VecGetArray(Y, &y));
858: for (i = 0; i < m; i++) y[i] = RDRadiationTemperature(rd, x[i].E);
859: PetscCall(VecRestoreArrayRead(X, (const PetscScalar **)&x));
860: PetscCall(VecRestoreArray(Y, &y));
862: PetscCall(VecView(Y, viewer));
863: PetscCall(VecDestroy(&Y));
864: PetscCall(DMDestroy(&da));
865: PetscFunctionReturn(PETSC_SUCCESS);
866: }
868: static PetscErrorCode RDTestDifferentiation(RD rd)
869: {
870: MPI_Comm comm;
871: RDNode n, nx;
872: PetscScalar epsilon;
874: PetscFunctionBeginUser;
875: PetscCall(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: PetscCall(PetscMemcpy(n1, n0, sizeof(n0)));
946: n1[i].E += epsilon;
947: fd[i].E = (RDDiffusion(rd, hx, n1, 1, 0) - a0) / epsilon;
948: PetscCall(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: PetscFunctionReturn(PETSC_SUCCESS);
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;
980: PetscFunctionBeginUser;
981: *inrd = 0;
982: PetscCall(PetscNew(&rd));
984: PetscOptionsBegin(comm, NULL, "Options for nonequilibrium radiation-diffusion with RD ionization", NULL);
985: {
986: rd->initial = 1;
987: PetscCall(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: PetscCall(PetscOptionsReal("-rd_unit_meter", "Length of 1 meter in nondimensional units", "", rd->unit.meter, &rd->unit.meter, 0));
1007: PetscCall(PetscOptionsReal("-rd_unit_kilogram", "Mass of 1 kilogram in nondimensional units", "", rd->unit.kilogram, &rd->unit.kilogram, 0));
1008: PetscCall(PetscOptionsReal("-rd_unit_second", "Time of a second in nondimensional units", "", rd->unit.second, &rd->unit.second, 0));
1009: PetscCall(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: PetscCall(PetscOptionsBool("-rd_monitor_residual", "Display residuals every time they are evaluated", "", rd->monitor_residual, &rd->monitor_residual, NULL));
1022: PetscCall(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: PetscCall(PetscOptionsEnum("-rd_quadrature", "Finite element quadrature", "", QuadratureTypes, (PetscEnum)rd->quadrature, (PetscEnum *)&rd->quadrature, NULL));
1026: }
1027: PetscCall(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: PetscCall(PetscOptionsEnum("-rd_leftbc", "Left boundary condition", "", BCTypes, (PetscEnum)rd->leftbc, (PetscEnum *)&rd->leftbc, NULL));
1057: PetscCall(PetscOptionsReal("-rd_E_applied", "Radiation flux at left end of domain", "", rd->Eapplied, &rd->Eapplied, NULL));
1058: PetscCall(PetscOptionsReal("-rd_beta", "Thermal exponent for photon absorption", "", rd->beta, &rd->beta, NULL));
1059: PetscCall(PetscOptionsReal("-rd_gamma", "Thermal exponent for diffusion coefficient", "", rd->gamma, &rd->gamma, NULL));
1060: PetscCall(PetscOptionsBool("-rd_view_draw", "Draw final solution", "", rd->view_draw, &rd->view_draw, NULL));
1061: PetscCall(PetscOptionsBool("-rd_endpoint", "Discretize using endpoints (like trapezoid rule) instead of midpoint", "", rd->endpoint, &rd->endpoint, NULL));
1062: PetscCall(PetscOptionsBool("-rd_bcmidpoint", "Impose the boundary condition at the midpoint (Theta) of the interval", "", rd->bcmidpoint, &rd->bcmidpoint, NULL));
1063: PetscCall(PetscOptionsBool("-rd_bclimit", "Limit diffusion coefficient in definition of Robin boundary condition", "", rd->bclimit, &rd->bclimit, NULL));
1064: PetscCall(PetscOptionsBool("-rd_test_diff", "Test differentiation in constitutive relations", "", rd->test_diff, &rd->test_diff, NULL));
1065: PetscCall(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: PetscCall(DMDACreate1d(comm, DM_BOUNDARY_NONE, 20, sizeof(RDNode) / sizeof(PetscScalar), 1, NULL, &rd->da));
1096: PetscCall(DMSetFromOptions(rd->da));
1097: PetscCall(DMSetUp(rd->da));
1098: PetscCall(DMDASetFieldName(rd->da, 0, "E"));
1099: PetscCall(DMDASetFieldName(rd->da, 1, "T"));
1100: PetscCall(DMDASetUniformCoordinates(rd->da, 0., 1., 0., 0., 0., 0.));
1102: *inrd = rd;
1103: PetscFunctionReturn(PETSC_SUCCESS);
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;
1116: PetscFunctionBeginUser;
1117: PetscCall(PetscInitialize(&argc, &argv, 0, NULL));
1118: PetscCall(RDCreate(PETSC_COMM_WORLD, &rd));
1119: PetscCall(DMCreateGlobalVector(rd->da, &X));
1120: PetscCall(DMSetMatType(rd->da, MATAIJ));
1121: PetscCall(DMCreateMatrix(rd->da, &B));
1122: PetscCall(RDInitialState(rd, X));
1124: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1125: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
1126: PetscCall(TSSetType(ts, TSTHETA));
1127: PetscCall(TSSetDM(ts, rd->da));
1128: switch (rd->discretization) {
1129: case DISCRETIZATION_FD:
1130: PetscCall(TSSetIFunction(ts, NULL, RDIFunction_FD, rd));
1131: if (rd->jacobian == JACOBIAN_ANALYTIC) PetscCall(TSSetIJacobian(ts, B, B, RDIJacobian_FD, rd));
1132: break;
1133: case DISCRETIZATION_FE:
1134: PetscCall(TSSetIFunction(ts, NULL, RDIFunction_FE, rd));
1135: if (rd->jacobian == JACOBIAN_ANALYTIC) PetscCall(TSSetIJacobian(ts, B, B, RDIJacobian_FE, rd));
1136: break;
1137: }
1138: PetscCall(TSSetMaxTime(ts, rd->final_time));
1139: PetscCall(TSSetTimeStep(ts, 1e-3));
1140: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1141: PetscCall(TSSetFromOptions(ts));
1143: A = B;
1144: PetscCall(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: PetscCall(SNESSetJacobian(snes, A, B, SNESComputeJacobianDefaultColor, 0));
1152: } break;
1153: case JACOBIAN_FD_FULL:
1154: PetscCall(SNESSetJacobian(snes, A, B, SNESComputeJacobianDefault, ts));
1155: break;
1156: }
1158: if (rd->test_diff) PetscCall(RDTestDifferentiation(rd));
1159: PetscCall(TSSolve(ts, X));
1160: PetscCall(TSGetSolveTime(ts, &ftime));
1161: PetscCall(TSGetStepNumber(ts, &steps));
1162: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Steps %" PetscInt_FMT " final time %g\n", steps, (double)ftime));
1163: if (rd->view_draw) PetscCall(RDView(rd, X, PETSC_VIEWER_DRAW_WORLD));
1164: if (rd->view_binary[0]) {
1165: PetscViewer viewer;
1166: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, rd->view_binary, FILE_MODE_WRITE, &viewer));
1167: PetscCall(RDView(rd, X, viewer));
1168: PetscCall(PetscViewerDestroy(&viewer));
1169: }
1170: PetscCall(VecDestroy(&X));
1171: PetscCall(MatDestroy(&B));
1172: PetscCall(RDDestroy(&rd));
1173: PetscCall(TSDestroy(&ts));
1174: PetscCall(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*/