Actual source code: ex10.c
1: static const char help[] = "1D nonequilibrium radiation diffusion with Saha ionization model.\n\n";
3: /*
4: This example implements the model described in
6: Rauenzahn, Mousseau, Knoll. "Temporal accuracy of the nonequilibrium radiation diffusion
7: equations employing a Saha ionization model" 2005.
9: The paper discusses three examples, the first two are nondimensional with a simple
10: ionization model. The third example is fully dimensional and uses the Saha ionization
11: model with realistic parameters.
12: */
14: #include <petscts.h>
15: #include <petscdm.h>
16: #include <petscdmda.h>
18: typedef enum {BC_DIRICHLET,BC_NEUMANN,BC_ROBIN} BCType;
19: static const char *const BCTypes[] = {"DIRICHLET","NEUMANN","ROBIN","BCType","BC_",0};
20: typedef enum {JACOBIAN_ANALYTIC,JACOBIAN_MATRIXFREE,JACOBIAN_FD_COLORING,JACOBIAN_FD_FULL} JacobianType;
21: static const char *const JacobianTypes[] = {"ANALYTIC","MATRIXFREE","FD_COLORING","FD_FULL","JacobianType","FD_",0};
22: typedef enum {DISCRETIZATION_FD,DISCRETIZATION_FE} DiscretizationType;
23: static const char *const DiscretizationTypes[] = {"FD","FE","DiscretizationType","DISCRETIZATION_",0};
24: typedef enum {QUADRATURE_GAUSS1,QUADRATURE_GAUSS2,QUADRATURE_GAUSS3,QUADRATURE_GAUSS4,QUADRATURE_LOBATTO2,QUADRATURE_LOBATTO3} QuadratureType;
25: static const char *const QuadratureTypes[] = {"GAUSS1","GAUSS2","GAUSS3","GAUSS4","LOBATTO2","LOBATTO3","QuadratureType","QUADRATURE_",0};
27: typedef struct {
28: PetscScalar E; /* radiation energy */
29: PetscScalar T; /* material temperature */
30: } RDNode;
32: typedef struct {
33: PetscReal meter,kilogram,second,Kelvin; /* Fundamental units */
34: PetscReal Joule,Watt; /* Derived units */
35: } RDUnit;
37: typedef struct _n_RD *RD;
39: struct _n_RD {
40: void (*MaterialEnergy)(RD,const RDNode*,PetscScalar*,RDNode*);
41: DM da;
42: PetscBool monitor_residual;
43: DiscretizationType discretization;
44: QuadratureType quadrature;
45: JacobianType jacobian;
46: PetscInt initial;
47: BCType leftbc;
48: PetscBool view_draw;
49: char view_binary[PETSC_MAX_PATH_LEN];
50: PetscBool test_diff;
51: PetscBool endpoint;
52: PetscBool bclimit;
53: PetscBool bcmidpoint;
54: RDUnit unit;
56: /* model constants, see Table 2 and RDCreate() */
57: PetscReal rho,K_R,K_p,I_H,m_p,m_e,h,k,c,sigma_b,beta,gamma;
59: /* Domain and boundary conditions */
60: PetscReal Eapplied; /* Radiation flux from the left */
61: PetscReal L; /* Length of domain */
62: PetscReal final_time;
63: };
65: static PetscErrorCode RDDestroy(RD *rd)
66: {
68: DMDestroy(&(*rd)->da);
69: PetscFree(*rd);
70: return 0;
71: }
73: /* The paper has a time derivative for material energy (Eq 2) which is a dependent variable (computable from temperature
74: * and density through an uninvertible relation). Computing this derivative is trivial for trapezoid rule (used in the
75: * paper), but does not generalize nicely to higher order integrators. Here we use the implicit form which provides
76: * time derivatives of the independent variables (radiation energy and temperature), so we must compute the time
77: * derivative of material energy ourselves (could be done using AD).
78: *
79: * There are multiple ionization models, this interface dispatches to the one currently in use.
80: */
81: static void RDMaterialEnergy(RD rd,const RDNode *n,PetscScalar *Em,RDNode *dEm) { rd->MaterialEnergy(rd,n,Em,dEm); }
83: /* Solves a quadratic equation while propagating tangents */
84: static void QuadraticSolve(PetscScalar a,PetscScalar a_t,PetscScalar b,PetscScalar b_t,PetscScalar c,PetscScalar c_t,PetscScalar *x,PetscScalar *x_t)
85: {
86: PetscScalar
87: disc = b*b - 4.*a*c,
88: disc_t = 2.*b*b_t - 4.*a_t*c - 4.*a*c_t,
89: num = -b + PetscSqrtScalar(disc), /* choose positive sign */
90: num_t = -b_t + 0.5/PetscSqrtScalar(disc)*disc_t,
91: den = 2.*a,
92: den_t = 2.*a_t;
93: *x = num/den;
94: *x_t = (num_t*den - num*den_t) / PetscSqr(den);
95: }
97: /* The primary model presented in the paper */
98: static void RDMaterialEnergy_Saha(RD rd,const RDNode *n,PetscScalar *inEm,RDNode *dEm)
99: {
100: PetscScalar Em,alpha,alpha_t,
101: T = n->T,
102: T_t = 1.,
103: chi = rd->I_H / (rd->k * T),
104: chi_t = -chi / T * T_t,
105: a = 1.,
106: a_t = 0,
107: 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 */
108: b_t = -b*chi_t + 1.5*b/chi*chi_t,
109: c = -b,
110: c_t = -b_t;
111: QuadraticSolve(a,a_t,b,b_t,c,c_t,&alpha,&alpha_t); /* Solve Eq 7 for alpha */
112: Em = rd->k * T / rd->m_p * (1.5*(1.+alpha) + alpha*chi); /* Eq 6 */
113: if (inEm) *inEm = Em;
114: if (dEm) {
115: dEm->E = 0;
116: dEm->T = Em / T * T_t + rd->k * T / rd->m_p * (1.5*alpha_t + alpha_t*chi + alpha*chi_t);
117: }
118: }
119: /* Reduced ionization model, Eq 30 */
120: static void RDMaterialEnergy_Reduced(RD rd,const RDNode *n,PetscScalar *Em,RDNode *dEm)
121: {
122: PetscScalar alpha,alpha_t,
123: T = n->T,
124: T_t = 1.,
125: chi = -0.3 / T,
126: chi_t = -chi / T * T_t,
127: a = 1.,
128: a_t = 0.,
129: b = PetscExpScalar(chi),
130: b_t = b*chi_t,
131: c = -b,
132: c_t = -b_t;
133: QuadraticSolve(a,a_t,b,b_t,c,c_t,&alpha,&alpha_t);
134: if (Em) *Em = (1.+alpha)*T + 0.3*alpha;
135: if (dEm) {
136: dEm->E = 0;
137: dEm->T = alpha_t*T + (1.+alpha)*T_t + 0.3*alpha_t;
138: }
139: }
141: /* Eq 5 */
142: static void RDSigma_R(RD rd,RDNode *n,PetscScalar *sigma_R,RDNode *dsigma_R)
143: {
144: *sigma_R = rd->K_R * rd->rho * PetscPowScalar(n->T,-rd->gamma);
145: dsigma_R->E = 0;
146: dsigma_R->T = -rd->gamma * (*sigma_R) / n->T;
147: }
149: /* Eq 4 */
150: static void RDDiffusionCoefficient(RD rd,PetscBool limit,RDNode *n,RDNode *nx,PetscScalar *D_R,RDNode *dD_R,RDNode *dxD_R)
151: {
152: PetscScalar sigma_R,denom;
153: RDNode dsigma_R,ddenom,dxdenom;
155: RDSigma_R(rd,n,&sigma_R,&dsigma_R);
156: denom = 3. * rd->rho * sigma_R + (int)limit * PetscAbsScalar(nx->E) / n->E;
157: ddenom.E = -(int)limit * PetscAbsScalar(nx->E) / PetscSqr(n->E);
158: ddenom.T = 3. * rd->rho * dsigma_R.T;
159: dxdenom.E = (int)limit * (PetscRealPart(nx->E)<0 ? -1. : 1.) / n->E;
160: dxdenom.T = 0;
161: *D_R = rd->c / denom;
162: if (dD_R) {
163: dD_R->E = -rd->c / PetscSqr(denom) * ddenom.E;
164: dD_R->T = -rd->c / PetscSqr(denom) * ddenom.T;
165: }
166: if (dxD_R) {
167: dxD_R->E = -rd->c / PetscSqr(denom) * dxdenom.E;
168: dxD_R->T = -rd->c / PetscSqr(denom) * dxdenom.T;
169: }
170: }
172: static PetscErrorCode RDStateView(RD rd,Vec X,Vec Xdot,Vec F)
173: {
175: DMDALocalInfo info;
176: PetscInt i;
177: const RDNode *x,*xdot,*f;
178: MPI_Comm comm;
181: PetscObjectGetComm((PetscObject)rd->da,&comm);
182: DMDAGetLocalInfo(rd->da,&info);
183: DMDAVecGetArrayRead(rd->da,X,(void*)&x);
184: DMDAVecGetArrayRead(rd->da,Xdot,(void*)&xdot);
185: DMDAVecGetArrayRead(rd->da,F,(void*)&f);
186: for (i=info.xs; i<info.xs+info.xm; i++) {
187: PetscSynchronizedPrintf(comm,"x[%D] (%10.2G,%10.2G) (%10.2G,%10.2G) (%10.2G,%10.2G)\n",i,PetscRealPart(x[i].E),PetscRealPart(x[i].T),
188: PetscRealPart(xdot[i].E),PetscRealPart(xdot[i].T), PetscRealPart(f[i].E),PetscRealPart(f[i].T));
189: }
190: DMDAVecRestoreArrayRead(rd->da,X,(void*)&x);
191: DMDAVecRestoreArrayRead(rd->da,Xdot,(void*)&xdot);
192: DMDAVecRestoreArrayRead(rd->da,F,(void*)&f);
193: PetscSynchronizedFlush(comm,PETSC_STDOUT);
194: return 0;
195: }
197: static PetscScalar RDRadiation(RD rd,const RDNode *n,RDNode *dn)
198: {
199: PetscScalar sigma_p = rd->K_p * rd->rho * PetscPowScalar(n->T,-rd->beta),
200: sigma_p_T = -rd->beta * sigma_p / n->T,
201: tmp = 4.* rd->sigma_b*PetscSqr(PetscSqr(n->T)) / rd->c - n->E,
202: tmp_E = -1.,
203: tmp_T = 4. * rd->sigma_b * 4 * n->T*(PetscSqr(n->T)) / rd->c,
204: rad = sigma_p * rd->c * rd->rho * tmp,
205: rad_E = sigma_p * rd->c * rd->rho * tmp_E,
206: rad_T = rd->c * rd->rho * (sigma_p_T * tmp + sigma_p * tmp_T);
207: if (dn) {
208: dn->E = rad_E;
209: dn->T = rad_T;
210: }
211: return rad;
212: }
214: static PetscScalar RDDiffusion(RD rd,PetscReal hx,const RDNode x[],PetscInt i,RDNode d[])
215: {
216: PetscReal ihx = 1./hx;
217: RDNode n_L,nx_L,n_R,nx_R,dD_L,dxD_L,dD_R,dxD_R,dfluxL[2],dfluxR[2];
218: PetscScalar D_L,D_R,fluxL,fluxR;
220: n_L.E = 0.5*(x[i-1].E + x[i].E);
221: n_L.T = 0.5*(x[i-1].T + x[i].T);
222: nx_L.E = (x[i].E - x[i-1].E)/hx;
223: nx_L.T = (x[i].T - x[i-1].T)/hx;
224: RDDiffusionCoefficient(rd,PETSC_TRUE,&n_L,&nx_L,&D_L,&dD_L,&dxD_L);
225: fluxL = D_L*nx_L.E;
226: dfluxL[0].E = -ihx*D_L + (0.5*dD_L.E - ihx*dxD_L.E)*nx_L.E;
227: dfluxL[1].E = +ihx*D_L + (0.5*dD_L.E + ihx*dxD_L.E)*nx_L.E;
228: dfluxL[0].T = (0.5*dD_L.T - ihx*dxD_L.T)*nx_L.E;
229: dfluxL[1].T = (0.5*dD_L.T + ihx*dxD_L.T)*nx_L.E;
231: n_R.E = 0.5*(x[i].E + x[i+1].E);
232: n_R.T = 0.5*(x[i].T + x[i+1].T);
233: nx_R.E = (x[i+1].E - x[i].E)/hx;
234: nx_R.T = (x[i+1].T - x[i].T)/hx;
235: RDDiffusionCoefficient(rd,PETSC_TRUE,&n_R,&nx_R,&D_R,&dD_R,&dxD_R);
236: fluxR = D_R*nx_R.E;
237: dfluxR[0].E = -ihx*D_R + (0.5*dD_R.E - ihx*dxD_R.E)*nx_R.E;
238: dfluxR[1].E = +ihx*D_R + (0.5*dD_R.E + ihx*dxD_R.E)*nx_R.E;
239: dfluxR[0].T = (0.5*dD_R.T - ihx*dxD_R.T)*nx_R.E;
240: dfluxR[1].T = (0.5*dD_R.T + ihx*dxD_R.T)*nx_R.E;
242: if (d) {
243: d[0].E = -ihx*dfluxL[0].E;
244: d[0].T = -ihx*dfluxL[0].T;
245: d[1].E = ihx*(dfluxR[0].E - dfluxL[1].E);
246: d[1].T = ihx*(dfluxR[0].T - dfluxL[1].T);
247: d[2].E = ihx*dfluxR[1].E;
248: d[2].T = ihx*dfluxR[1].T;
249: }
250: return ihx*(fluxR - fluxL);
251: }
253: 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)
254: {
255: PetscBool istheta;
258: DMGetLocalVector(rd->da,X0loc);
259: DMGetLocalVector(rd->da,Xloc);
260: DMGetLocalVector(rd->da,Xloc_t);
262: DMGlobalToLocalBegin(rd->da,X,INSERT_VALUES,*Xloc);
263: DMGlobalToLocalEnd(rd->da,X,INSERT_VALUES,*Xloc);
264: DMGlobalToLocalBegin(rd->da,Xdot,INSERT_VALUES,*Xloc_t);
265: DMGlobalToLocalEnd(rd->da,Xdot,INSERT_VALUES,*Xloc_t);
267: /*
268: The following is a hack to subvert TSTHETA which is like an implicit midpoint method to behave more like a trapezoid
269: rule. These methods have equivalent linear stability, but the nonlinear stability is somewhat different. The
270: radiation system is inconvenient to write in explicit form because the ionization model is "on the left".
271: */
272: PetscObjectTypeCompare((PetscObject)ts,TSTHETA,&istheta);
273: if (istheta && rd->endpoint) {
274: TSThetaGetTheta(ts,Theta);
275: } else *Theta = 1.;
277: TSGetTimeStep(ts,dt);
278: VecWAXPY(*X0loc,-(*Theta)*(*dt),*Xloc_t,*Xloc); /* back out the value at the start of this step */
279: if (rd->endpoint) {
280: VecWAXPY(*Xloc,*dt,*Xloc_t,*X0loc); /* move the abscissa to the end of the step */
281: }
283: DMDAVecGetArray(rd->da,*X0loc,x0);
284: DMDAVecGetArray(rd->da,*Xloc,x);
285: DMDAVecGetArray(rd->da,*Xloc_t,xdot);
286: return 0;
287: }
289: static PetscErrorCode RDRestoreLocalArrays(RD rd,Vec *X0loc,RDNode **x0,Vec *Xloc,RDNode **x,Vec *Xloc_t,RDNode **xdot)
290: {
292: DMDAVecRestoreArray(rd->da,*X0loc,x0);
293: DMDAVecRestoreArray(rd->da,*Xloc,x);
294: DMDAVecRestoreArray(rd->da,*Xloc_t,xdot);
295: DMRestoreLocalVector(rd->da,X0loc);
296: DMRestoreLocalVector(rd->da,Xloc);
297: DMRestoreLocalVector(rd->da,Xloc_t);
298: return 0;
299: }
301: static PetscErrorCode PETSC_UNUSED RDCheckDomain_Private(RD rd,TS ts,Vec X,PetscBool *in)
302: {
303: PetscInt minloc;
304: PetscReal min;
307: VecMin(X,&minloc,&min);
308: if (min < 0) {
309: SNES snes;
310: *in = PETSC_FALSE;
311: TSGetSNES(ts,&snes);
312: SNESSetFunctionDomainError(snes);
313: PetscInfo(ts,"Domain violation at %D field %D value %g\n",minloc/2,minloc%2,(double)min);
314: } else *in = PETSC_TRUE;
315: return 0;
316: }
318: /* Energy and temperature must remain positive */
319: #define RDCheckDomain(rd,ts,X) do { \
320: PetscBool _in; \
321: RDCheckDomain_Private(rd,ts,X,&_in); \
322: if (!_in) return 0; \
323: } while (0)
325: static PetscErrorCode RDIFunction_FD(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
326: {
327: RD rd = (RD)ctx;
328: RDNode *x,*x0,*xdot,*f;
329: Vec X0loc,Xloc,Xloc_t;
330: PetscReal hx,Theta,dt;
331: DMDALocalInfo info;
332: PetscInt i;
335: RDGetLocalArrays(rd,ts,X,Xdot,&Theta,&dt,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
336: DMDAVecGetArray(rd->da,F,&f);
337: DMDAGetLocalInfo(rd->da,&info);
338: VecZeroEntries(F);
340: hx = rd->L / (info.mx-1);
342: for (i=info.xs; i<info.xs+info.xm; i++) {
343: PetscReal rho = rd->rho;
344: PetscScalar Em_t,rad;
346: rad = (1.-Theta)*RDRadiation(rd,&x0[i],0) + Theta*RDRadiation(rd,&x[i],0);
347: if (rd->endpoint) {
348: PetscScalar Em0,Em1;
349: RDMaterialEnergy(rd,&x0[i],&Em0,NULL);
350: RDMaterialEnergy(rd,&x[i],&Em1,NULL);
351: Em_t = (Em1 - Em0) / dt;
352: } else {
353: RDNode dEm;
354: RDMaterialEnergy(rd,&x[i],NULL,&dEm);
355: Em_t = dEm.E * xdot[i].E + dEm.T * xdot[i].T;
356: }
357: /* Residuals are multiplied by the volume element (hx). */
358: /* The temperature equation does not have boundary conditions */
359: f[i].T = hx*(rho*Em_t + rad);
361: if (i == 0) { /* Left boundary condition */
362: PetscScalar D_R,bcTheta = rd->bcmidpoint ? Theta : 1.;
363: RDNode n, nx;
365: n.E = (1.-bcTheta)*x0[0].E + bcTheta*x[0].E;
366: n.T = (1.-bcTheta)*x0[0].T + bcTheta*x[0].T;
367: nx.E = ((1.-bcTheta)*(x0[1].E-x0[0].E) + bcTheta*(x[1].E-x[0].E))/hx;
368: nx.T = ((1.-bcTheta)*(x0[1].T-x0[0].T) + bcTheta*(x[1].T-x[0].T))/hx;
369: switch (rd->leftbc) {
370: case BC_ROBIN:
371: RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D_R,0,0);
372: f[0].E = hx*(n.E - 2. * D_R * nx.E - rd->Eapplied);
373: break;
374: case BC_NEUMANN:
375: f[0].E = x[1].E - x[0].E;
376: break;
377: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D",rd->initial);
378: }
379: } else if (i == info.mx-1) { /* Right boundary */
380: f[i].E = x[i].E - x[i-1].E; /* Homogeneous Neumann */
381: } else {
382: PetscScalar diff = (1.-Theta)*RDDiffusion(rd,hx,x0,i,0) + Theta*RDDiffusion(rd,hx,x,i,0);
383: f[i].E = hx*(xdot[i].E - diff - rad);
384: }
385: }
386: RDRestoreLocalArrays(rd,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
387: DMDAVecRestoreArray(rd->da,F,&f);
388: if (rd->monitor_residual) RDStateView(rd,X,Xdot,F);
389: return 0;
390: }
392: static PetscErrorCode RDIJacobian_FD(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
393: {
394: RD rd = (RD)ctx;
395: RDNode *x,*x0,*xdot;
396: Vec X0loc,Xloc,Xloc_t;
397: PetscReal hx,Theta,dt;
398: DMDALocalInfo info;
399: PetscInt i;
402: RDGetLocalArrays(rd,ts,X,Xdot,&Theta,&dt,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
403: DMDAGetLocalInfo(rd->da,&info);
404: hx = rd->L / (info.mx-1);
405: MatZeroEntries(B);
407: for (i=info.xs; i<info.xs+info.xm; i++) {
408: PetscInt col[3];
409: PetscReal rho = rd->rho;
410: PetscScalar /*Em_t,rad,*/ K[2][6];
411: RDNode dEm_t,drad;
413: /*rad = (1.-Theta)* */ RDRadiation(rd,&x0[i],0); /* + Theta* */ RDRadiation(rd,&x[i],&drad);
415: if (rd->endpoint) {
416: PetscScalar Em0,Em1;
417: RDNode dEm1;
418: RDMaterialEnergy(rd,&x0[i],&Em0,NULL);
419: RDMaterialEnergy(rd,&x[i],&Em1,&dEm1);
420: /*Em_t = (Em1 - Em0) / (Theta*dt);*/
421: dEm_t.E = dEm1.E / (Theta*dt);
422: dEm_t.T = dEm1.T / (Theta*dt);
423: } else {
424: const PetscScalar epsilon = x[i].T * PETSC_SQRT_MACHINE_EPSILON;
425: RDNode n1;
426: RDNode dEm,dEm1;
427: PetscScalar Em_TT;
429: n1.E = x[i].E;
430: n1.T = x[i].T+epsilon;
431: RDMaterialEnergy(rd,&x[i],NULL,&dEm);
432: RDMaterialEnergy(rd,&n1,NULL,&dEm1);
433: /* The Jacobian needs another derivative. We finite difference here instead of
434: * propagating second derivatives through the ionization model. */
435: Em_TT = (dEm1.T - dEm.T) / epsilon;
436: /*Em_t = dEm.E * xdot[i].E + dEm.T * xdot[i].T;*/
437: dEm_t.E = dEm.E * a;
438: dEm_t.T = dEm.T * a + Em_TT * xdot[i].T;
439: }
441: PetscMemzero(K,sizeof(K));
442: /* Residuals are multiplied by the volume element (hx). */
443: if (i == 0) {
444: PetscScalar D,bcTheta = rd->bcmidpoint ? Theta : 1.;
445: RDNode n, nx;
446: RDNode dD,dxD;
448: n.E = (1.-bcTheta)*x0[0].E + bcTheta*x[0].E;
449: n.T = (1.-bcTheta)*x0[0].T + bcTheta*x[0].T;
450: nx.E = ((1.-bcTheta)*(x0[1].E-x0[0].E) + bcTheta*(x[1].E-x[0].E))/hx;
451: nx.T = ((1.-bcTheta)*(x0[1].T-x0[0].T) + bcTheta*(x[1].T-x[0].T))/hx;
452: switch (rd->leftbc) {
453: case BC_ROBIN:
454: RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D,&dD,&dxD);
455: K[0][1*2+0] = (bcTheta/Theta)*hx*(1. -2.*D*(-1./hx) - 2.*nx.E*dD.E + 2.*nx.E*dxD.E/hx);
456: K[0][1*2+1] = (bcTheta/Theta)*hx*(-2.*nx.E*dD.T);
457: K[0][2*2+0] = (bcTheta/Theta)*hx*(-2.*D*(1./hx) - 2.*nx.E*dD.E - 2.*nx.E*dxD.E/hx);
458: break;
459: case BC_NEUMANN:
460: K[0][1*2+0] = -1./Theta;
461: K[0][2*2+0] = 1./Theta;
462: break;
463: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D",rd->initial);
464: }
465: } else if (i == info.mx-1) {
466: K[0][0*2+0] = -1./Theta;
467: K[0][1*2+0] = 1./Theta;
468: } else {
469: /*PetscScalar diff;*/
470: RDNode ddiff[3];
471: /*diff = (1.-Theta)*RDDiffusion(rd,hx,x0,i,0) + Theta* */ RDDiffusion(rd,hx,x,i,ddiff);
472: K[0][0*2+0] = -hx*ddiff[0].E;
473: K[0][0*2+1] = -hx*ddiff[0].T;
474: K[0][1*2+0] = hx*(a - ddiff[1].E - drad.E);
475: K[0][1*2+1] = hx*(-ddiff[1].T - drad.T);
476: K[0][2*2+0] = -hx*ddiff[2].E;
477: K[0][2*2+1] = -hx*ddiff[2].T;
478: }
480: K[1][1*2+0] = hx*(rho*dEm_t.E + drad.E);
481: K[1][1*2+1] = hx*(rho*dEm_t.T + drad.T);
483: col[0] = i-1;
484: col[1] = i;
485: col[2] = i+1<info.mx ? i+1 : -1;
486: MatSetValuesBlocked(B,1,&i,3,col,&K[0][0],INSERT_VALUES);
487: }
488: RDRestoreLocalArrays(rd,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
489: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
490: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
491: if (A != B) {
492: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
493: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
494: }
495: return 0;
496: }
498: /* Evaluate interpolants and derivatives at a select quadrature point */
499: static void RDEvaluate(PetscReal interp[][2],PetscReal deriv[][2],PetscInt q,const RDNode x[],PetscInt i,RDNode *n,RDNode *nx)
500: {
501: PetscInt j;
502: n->E = 0; n->T = 0; nx->E = 0; nx->T = 0;
503: for (j=0; j<2; j++) {
504: n->E += interp[q][j] * x[i+j].E;
505: n->T += interp[q][j] * x[i+j].T;
506: nx->E += deriv[q][j] * x[i+j].E;
507: nx->T += deriv[q][j] * x[i+j].T;
508: }
509: }
511: /*
512: Various quadrature rules. The nonlinear terms are non-polynomial so no standard quadrature will be exact.
513: */
514: static PetscErrorCode RDGetQuadrature(RD rd,PetscReal hx,PetscInt *nq,PetscReal weight[],PetscReal interp[][2],PetscReal deriv[][2])
515: {
516: PetscInt q,j;
517: const PetscReal *refweight,(*refinterp)[2],(*refderiv)[2];
520: switch (rd->quadrature) {
521: case QUADRATURE_GAUSS1: {
522: static const PetscReal ww[1] = {1.},ii[1][2] = {{0.5,0.5}},dd[1][2] = {{-1.,1.}};
523: *nq = 1; refweight = ww; refinterp = ii; refderiv = dd;
524: } break;
525: case QUADRATURE_GAUSS2: {
526: static const PetscReal ii[2][2] = {{0.78867513459481287,0.21132486540518713},{0.21132486540518713,0.78867513459481287}},dd[2][2] = {{-1.,1.},{-1.,1.}},ww[2] = {0.5,0.5};
527: *nq = 2; refweight = ww; refinterp = ii; refderiv = dd;
528: } break;
529: case QUADRATURE_GAUSS3: {
530: static const PetscReal ii[3][2] = {{0.8872983346207417,0.1127016653792583},{0.5,0.5},{0.1127016653792583,0.8872983346207417}},
531: dd[3][2] = {{-1,1},{-1,1},{-1,1}},ww[3] = {5./18,8./18,5./18};
532: *nq = 3; refweight = ww; refinterp = ii; refderiv = dd;
533: } break;
534: case QUADRATURE_GAUSS4: {
535: static const PetscReal ii[][2] = {{0.93056815579702623,0.069431844202973658},
536: {0.66999052179242813,0.33000947820757187},
537: {0.33000947820757187,0.66999052179242813},
538: {0.069431844202973658,0.93056815579702623}},
539: dd[][2] = {{-1,1},{-1,1},{-1,1},{-1,1}},ww[] = {0.17392742256872692,0.3260725774312731,0.3260725774312731,0.17392742256872692};
541: *nq = 4; refweight = ww; refinterp = ii; refderiv = dd;
542: } break;
543: case QUADRATURE_LOBATTO2: {
544: static const PetscReal ii[2][2] = {{1.,0.},{0.,1.}},dd[2][2] = {{-1.,1.},{-1.,1.}},ww[2] = {0.5,0.5};
545: *nq = 2; refweight = ww; refinterp = ii; refderiv = dd;
546: } break;
547: case QUADRATURE_LOBATTO3: {
548: static const PetscReal ii[3][2] = {{1,0},{0.5,0.5},{0,1}},dd[3][2] = {{-1,1},{-1,1},{-1,1}},ww[3] = {1./6,4./6,1./6};
549: *nq = 3; refweight = ww; refinterp = ii; refderiv = dd;
550: } break;
551: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown quadrature %d",(int)rd->quadrature);
552: }
554: for (q=0; q<*nq; q++) {
555: weight[q] = refweight[q] * hx;
556: for (j=0; j<2; j++) {
557: interp[q][j] = refinterp[q][j];
558: deriv[q][j] = refderiv[q][j] / hx;
559: }
560: }
561: return 0;
562: }
564: /*
565: Finite element version
566: */
567: static PetscErrorCode RDIFunction_FE(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
568: {
569: RD rd = (RD)ctx;
570: RDNode *x,*x0,*xdot,*f;
571: Vec X0loc,Xloc,Xloc_t,Floc;
572: PetscReal hx,Theta,dt,weight[5],interp[5][2],deriv[5][2];
573: DMDALocalInfo info;
574: PetscInt i,j,q,nq;
577: RDGetLocalArrays(rd,ts,X,Xdot,&Theta,&dt,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
579: DMGetLocalVector(rd->da,&Floc);
580: VecZeroEntries(Floc);
581: DMDAVecGetArray(rd->da,Floc,&f);
582: DMDAGetLocalInfo(rd->da,&info);
584: /* Set up shape functions and quadrature for elements (assumes a uniform grid) */
585: hx = rd->L / (info.mx-1);
586: RDGetQuadrature(rd,hx,&nq,weight,interp,deriv);
588: for (i=info.xs; i<PetscMin(info.xs+info.xm,info.mx-1); i++) {
589: for (q=0; q<nq; q++) {
590: PetscReal rho = rd->rho;
591: PetscScalar Em_t,rad,D_R,D0_R;
592: RDNode n,n0,nx,n0x,nt,ntx;
593: RDEvaluate(interp,deriv,q,x,i,&n,&nx);
594: RDEvaluate(interp,deriv,q,x0,i,&n0,&n0x);
595: RDEvaluate(interp,deriv,q,xdot,i,&nt,&ntx);
597: rad = (1.-Theta)*RDRadiation(rd,&n0,0) + Theta*RDRadiation(rd,&n,0);
598: if (rd->endpoint) {
599: PetscScalar Em0,Em1;
600: RDMaterialEnergy(rd,&n0,&Em0,NULL);
601: RDMaterialEnergy(rd,&n,&Em1,NULL);
602: Em_t = (Em1 - Em0) / dt;
603: } else {
604: RDNode dEm;
605: RDMaterialEnergy(rd,&n,NULL,&dEm);
606: Em_t = dEm.E * nt.E + dEm.T * nt.T;
607: }
608: RDDiffusionCoefficient(rd,PETSC_TRUE,&n0,&n0x,&D0_R,0,0);
609: RDDiffusionCoefficient(rd,PETSC_TRUE,&n,&nx,&D_R,0,0);
610: for (j=0; j<2; j++) {
611: f[i+j].E += (deriv[q][j] * weight[q] * ((1.-Theta)*D0_R*n0x.E + Theta*D_R*nx.E)
612: + interp[q][j] * weight[q] * (nt.E - rad));
613: f[i+j].T += interp[q][j] * weight[q] * (rho * Em_t + rad);
614: }
615: }
616: }
617: if (info.xs == 0) {
618: switch (rd->leftbc) {
619: case BC_ROBIN: {
620: PetscScalar D_R,D_R_bc;
621: PetscReal ratio,bcTheta = rd->bcmidpoint ? Theta : 1.;
622: RDNode n, nx;
624: n.E = (1-bcTheta)*x0[0].E + bcTheta*x[0].E;
625: n.T = (1-bcTheta)*x0[0].T + bcTheta*x[0].T;
626: nx.E = (x[1].E-x[0].E)/hx;
627: nx.T = (x[1].T-x[0].T)/hx;
628: RDDiffusionCoefficient(rd,PETSC_TRUE,&n,&nx,&D_R,0,0);
629: RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D_R_bc,0,0);
630: ratio = PetscRealPart(D_R/D_R_bc);
633: f[0].E += -ratio*0.5*(rd->Eapplied - n.E);
634: } break;
635: case BC_NEUMANN:
636: /* homogeneous Neumann is the natural condition */
637: break;
638: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D",rd->initial);
639: }
640: }
642: RDRestoreLocalArrays(rd,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
643: DMDAVecRestoreArray(rd->da,Floc,&f);
644: VecZeroEntries(F);
645: DMLocalToGlobalBegin(rd->da,Floc,ADD_VALUES,F);
646: DMLocalToGlobalEnd(rd->da,Floc,ADD_VALUES,F);
647: DMRestoreLocalVector(rd->da,&Floc);
649: if (rd->monitor_residual) RDStateView(rd,X,Xdot,F);
650: return 0;
651: }
653: static PetscErrorCode RDIJacobian_FE(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
654: {
655: RD rd = (RD)ctx;
656: RDNode *x,*x0,*xdot;
657: Vec X0loc,Xloc,Xloc_t;
658: PetscReal hx,Theta,dt,weight[5],interp[5][2],deriv[5][2];
659: DMDALocalInfo info;
660: PetscInt i,j,k,q,nq;
661: PetscScalar K[4][4];
664: RDGetLocalArrays(rd,ts,X,Xdot,&Theta,&dt,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
665: DMDAGetLocalInfo(rd->da,&info);
666: hx = rd->L / (info.mx-1);
667: RDGetQuadrature(rd,hx,&nq,weight,interp,deriv);
668: MatZeroEntries(B);
669: for (i=info.xs; i<PetscMin(info.xs+info.xm,info.mx-1); i++) {
670: PetscInt rc[2];
672: rc[0] = i; rc[1] = i+1;
673: PetscMemzero(K,sizeof(K));
674: for (q=0; q<nq; q++) {
675: PetscScalar D_R;
676: PETSC_UNUSED PetscScalar rad;
677: RDNode n,nx,nt,ntx,drad,dD_R,dxD_R,dEm;
678: RDEvaluate(interp,deriv,q,x,i,&n,&nx);
679: RDEvaluate(interp,deriv,q,xdot,i,&nt,&ntx);
680: rad = RDRadiation(rd,&n,&drad);
681: RDDiffusionCoefficient(rd,PETSC_TRUE,&n,&nx,&D_R,&dD_R,&dxD_R);
682: RDMaterialEnergy(rd,&n,NULL,&dEm);
683: for (j=0; j<2; j++) {
684: for (k=0; k<2; k++) {
685: K[j*2+0][k*2+0] += (+interp[q][j] * weight[q] * (a - drad.E) * interp[q][k]
686: + deriv[q][j] * weight[q] * ((D_R + dxD_R.E * nx.E) * deriv[q][k] + dD_R.E * nx.E * interp[q][k]));
687: K[j*2+0][k*2+1] += (+interp[q][j] * weight[q] * (-drad.T * interp[q][k])
688: + deriv[q][j] * weight[q] * (dxD_R.T * deriv[q][k] + dD_R.T * interp[q][k]) * nx.E);
689: K[j*2+1][k*2+0] += interp[q][j] * weight[q] * drad.E * interp[q][k];
690: K[j*2+1][k*2+1] += interp[q][j] * weight[q] * (a * rd->rho * dEm.T + drad.T) * interp[q][k];
691: }
692: }
693: }
694: MatSetValuesBlocked(B,2,rc,2,rc,&K[0][0],ADD_VALUES);
695: }
696: if (info.xs == 0) {
697: switch (rd->leftbc) {
698: case BC_ROBIN: {
699: PetscScalar D_R,D_R_bc;
700: PetscReal ratio;
701: RDNode n, nx;
703: n.E = (1-Theta)*x0[0].E + Theta*x[0].E;
704: n.T = (1-Theta)*x0[0].T + Theta*x[0].T;
705: nx.E = (x[1].E-x[0].E)/hx;
706: nx.T = (x[1].T-x[0].T)/hx;
707: RDDiffusionCoefficient(rd,PETSC_TRUE,&n,&nx,&D_R,0,0);
708: RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D_R_bc,0,0);
709: ratio = PetscRealPart(D_R/D_R_bc);
710: MatSetValue(B,0,0,ratio*0.5,ADD_VALUES);
711: } break;
712: case BC_NEUMANN:
713: /* homogeneous Neumann is the natural condition */
714: break;
715: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D",rd->initial);
716: }
717: }
719: RDRestoreLocalArrays(rd,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
720: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
721: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
722: if (A != B) {
723: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
724: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
725: }
726: return 0;
727: }
729: /* Temperature that is in equilibrium with the radiation density */
730: static PetscScalar RDRadiationTemperature(RD rd,PetscScalar E) { return PetscPowScalar(E*rd->c/(4.*rd->sigma_b),0.25); }
732: static PetscErrorCode RDInitialState(RD rd,Vec X)
733: {
734: DMDALocalInfo info;
735: PetscInt i;
736: RDNode *x;
739: DMDAGetLocalInfo(rd->da,&info);
740: DMDAVecGetArray(rd->da,X,&x);
741: for (i=info.xs; i<info.xs+info.xm; i++) {
742: PetscReal coord = i*rd->L/(info.mx-1);
743: switch (rd->initial) {
744: case 1:
745: x[i].E = 0.001;
746: x[i].T = RDRadiationTemperature(rd,x[i].E);
747: break;
748: case 2:
749: x[i].E = 0.001 + 100.*PetscExpReal(-PetscSqr(coord/0.1));
750: x[i].T = RDRadiationTemperature(rd,x[i].E);
751: break;
752: case 3:
753: x[i].E = 7.56e-2 * rd->unit.Joule / PetscPowScalarInt(rd->unit.meter,3);
754: x[i].T = RDRadiationTemperature(rd,x[i].E);
755: break;
756: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No initial state %D",rd->initial);
757: }
758: }
759: DMDAVecRestoreArray(rd->da,X,&x);
760: return 0;
761: }
763: static PetscErrorCode RDView(RD rd,Vec X,PetscViewer viewer)
764: {
765: Vec Y;
766: const RDNode *x;
767: PetscScalar *y;
768: PetscInt i,m,M;
769: const PetscInt *lx;
770: DM da;
771: MPI_Comm comm;
774: /*
775: Create a DMDA (one dof per node, zero stencil width, same layout) to hold Trad
776: (radiation temperature). It is not necessary to create a DMDA for this, but this way
777: output and visualization will have meaningful variable names and correct scales.
778: */
779: DMDAGetInfo(rd->da,0, &M,0,0, 0,0,0, 0,0,0,0,0,0);
780: DMDAGetOwnershipRanges(rd->da,&lx,0,0);
781: PetscObjectGetComm((PetscObject)rd->da,&comm);
782: DMDACreate1d(comm,DM_BOUNDARY_NONE,M,1,0,lx,&da);
783: DMSetFromOptions(da);
784: DMSetUp(da);
785: DMDASetUniformCoordinates(da,0.,rd->L,0.,0.,0.,0.);
786: DMDASetFieldName(da,0,"T_rad");
787: DMCreateGlobalVector(da,&Y);
789: /* Compute the radiation temperature from the solution at each node */
790: VecGetLocalSize(Y,&m);
791: VecGetArrayRead(X,(const PetscScalar **)&x);
792: VecGetArray(Y,&y);
793: for (i=0; i<m; i++) y[i] = RDRadiationTemperature(rd,x[i].E);
794: VecRestoreArrayRead(X,(const PetscScalar**)&x);
795: VecRestoreArray(Y,&y);
797: VecView(Y,viewer);
798: VecDestroy(&Y);
799: DMDestroy(&da);
800: return 0;
801: }
803: static PetscErrorCode RDTestDifferentiation(RD rd)
804: {
805: MPI_Comm comm;
807: RDNode n,nx;
808: PetscScalar epsilon;
811: PetscObjectGetComm((PetscObject)rd->da,&comm);
812: epsilon = 1e-8;
813: {
814: RDNode dEm,fdEm;
815: PetscScalar T0 = 1000.,T1 = T0*(1.+epsilon),Em0,Em1;
816: n.E = 1.;
817: n.T = T0;
818: rd->MaterialEnergy(rd,&n,&Em0,&dEm);
819: n.E = 1.+epsilon;
820: n.T = T0;
821: rd->MaterialEnergy(rd,&n,&Em1,0);
822: fdEm.E = (Em1-Em0)/epsilon;
823: n.E = 1.;
824: n.T = T1;
825: rd->MaterialEnergy(rd,&n,&Em1,0);
826: fdEm.T = (Em1-Em0)/(T0*epsilon);
827: PetscPrintf(comm,"dEm {%g,%g}, fdEm {%g,%g}, diff {%g,%g}\n",(double)PetscRealPart(dEm.E),(double)PetscRealPart(dEm.T),
828: (double)PetscRealPart(fdEm.E),(double)PetscRealPart(fdEm.T),(double)PetscRealPart(dEm.E-fdEm.E),(double)PetscRealPart(dEm.T-fdEm.T));
829: }
830: {
831: PetscScalar D0,D;
832: RDNode dD,dxD,fdD,fdxD;
833: n.E = 1.; n.T = 1.; nx.E = 1.; n.T = 1.;
834: RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D0,&dD,&dxD);
835: n.E = 1.+epsilon; n.T = 1.; nx.E = 1.; n.T = 1.;
836: RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D,0,0); fdD.E = (D-D0)/epsilon;
837: n.E = 1; n.T = 1.+epsilon; nx.E = 1.; n.T = 1.;
838: RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D,0,0); fdD.T = (D-D0)/epsilon;
839: n.E = 1; n.T = 1.; nx.E = 1.+epsilon; n.T = 1.;
840: RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D,0,0); fdxD.E = (D-D0)/epsilon;
841: n.E = 1; n.T = 1.; nx.E = 1.; n.T = 1.+epsilon;
842: RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D,0,0); fdxD.T = (D-D0)/epsilon;
843: PetscPrintf(comm,"dD {%g,%g}, fdD {%g,%g}, diff {%g,%g}\n",(double)PetscRealPart(dD.E),(double)PetscRealPart(dD.T),
844: (double)PetscRealPart(fdD.E),(double)PetscRealPart(fdD.T),(double)PetscRealPart(dD.E-fdD.E),(double)PetscRealPart(dD.T-fdD.T));
845: PetscPrintf(comm,"dxD {%g,%g}, fdxD {%g,%g}, diffx {%g,%g}\n",(double)PetscRealPart(dxD.E),(double)PetscRealPart(dxD.T),
846: (double)PetscRealPart(fdxD.E),(double)PetscRealPart(fdxD.T),(double)PetscRealPart(dxD.E-fdxD.E),(double)PetscRealPart(dxD.T-fdxD.T));
847: }
848: {
849: PetscInt i;
850: PetscReal hx = 1.;
851: PetscScalar a0;
852: RDNode n0[3],n1[3],d[3],fd[3];
854: n0[0].E = 1.;
855: n0[0].T = 1.;
856: n0[1].E = 5.;
857: n0[1].T = 3.;
858: n0[2].E = 4.;
859: n0[2].T = 2.;
860: a0 = RDDiffusion(rd,hx,n0,1,d);
861: for (i=0; i<3; i++) {
862: PetscMemcpy(n1,n0,sizeof(n0)); n1[i].E += epsilon;
863: fd[i].E = (RDDiffusion(rd,hx,n1,1,0)-a0)/epsilon;
864: PetscMemcpy(n1,n0,sizeof(n0)); n1[i].T += epsilon;
865: fd[i].T = (RDDiffusion(rd,hx,n1,1,0)-a0)/epsilon;
866: PetscPrintf(comm,"ddiff[%D] {%g,%g}, fd {%g %g}, diff {%g,%g}\n",i,(double)PetscRealPart(d[i].E),(double)PetscRealPart(d[i].T),
867: (double)PetscRealPart(fd[i].E),(double)PetscRealPart(fd[i].T),(double)PetscRealPart(d[i].E-fd[i].E),(double)PetscRealPart(d[i].T-fd[i].T));
868: }
869: }
870: {
871: PetscScalar rad0,rad;
872: RDNode drad,fdrad;
873: n.E = 1.; n.T = 1.;
874: rad0 = RDRadiation(rd,&n,&drad);
875: n.E = 1.+epsilon; n.T = 1.;
876: rad = RDRadiation(rd,&n,0); fdrad.E = (rad-rad0)/epsilon;
877: n.E = 1.; n.T = 1.+epsilon;
878: rad = RDRadiation(rd,&n,0); fdrad.T = (rad-rad0)/epsilon;
879: PetscPrintf(comm,"drad {%g,%g}, fdrad {%g,%g}, diff {%g,%g}\n",(double)PetscRealPart(drad.E),(double)PetscRealPart(drad.T),
880: (double)PetscRealPart(fdrad.E),(double)PetscRealPart(fdrad.T),(double)PetscRealPart(drad.E-drad.E),(double)PetscRealPart(drad.T-fdrad.T));
881: }
882: return 0;
883: }
885: static PetscErrorCode RDCreate(MPI_Comm comm,RD *inrd)
886: {
888: RD rd;
889: PetscReal meter=0,kilogram=0,second=0,Kelvin=0,Joule=0,Watt=0;
892: *inrd = 0;
893: PetscNew(&rd);
895: PetscOptionsBegin(comm,NULL,"Options for nonequilibrium radiation-diffusion with RD ionization",NULL);
896: {
897: rd->initial = 1;
898: PetscOptionsInt("-rd_initial","Initial condition (1=Marshak, 2=Blast, 3=Marshak+)","",rd->initial,&rd->initial,0);
899: switch (rd->initial) {
900: case 1:
901: case 2:
902: rd->unit.kilogram = 1.;
903: rd->unit.meter = 1.;
904: rd->unit.second = 1.;
905: rd->unit.Kelvin = 1.;
906: break;
907: case 3:
908: rd->unit.kilogram = 1.e12;
909: rd->unit.meter = 1.;
910: rd->unit.second = 1.e9;
911: rd->unit.Kelvin = 1.;
912: break;
913: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown initial condition %d",rd->initial);
914: }
915: /* Fundamental units */
916: PetscOptionsReal("-rd_unit_meter","Length of 1 meter in nondimensional units","",rd->unit.meter,&rd->unit.meter,0);
917: PetscOptionsReal("-rd_unit_kilogram","Mass of 1 kilogram in nondimensional units","",rd->unit.kilogram,&rd->unit.kilogram,0);
918: PetscOptionsReal("-rd_unit_second","Time of a second in nondimensional units","",rd->unit.second,&rd->unit.second,0);
919: PetscOptionsReal("-rd_unit_Kelvin","Temperature of a Kelvin in nondimensional units","",rd->unit.Kelvin,&rd->unit.Kelvin,0);
920: /* Derived units */
921: rd->unit.Joule = rd->unit.kilogram*PetscSqr(rd->unit.meter/rd->unit.second);
922: rd->unit.Watt = rd->unit.Joule/rd->unit.second;
923: /* Local aliases */
924: meter = rd->unit.meter;
925: kilogram = rd->unit.kilogram;
926: second = rd->unit.second;
927: Kelvin = rd->unit.Kelvin;
928: Joule = rd->unit.Joule;
929: Watt = rd->unit.Watt;
931: PetscOptionsBool("-rd_monitor_residual","Display residuals every time they are evaluated","",rd->monitor_residual,&rd->monitor_residual,NULL);
932: PetscOptionsEnum("-rd_discretization","Discretization type","",DiscretizationTypes,(PetscEnum)rd->discretization,(PetscEnum*)&rd->discretization,NULL);
933: if (rd->discretization == DISCRETIZATION_FE) {
934: rd->quadrature = QUADRATURE_GAUSS2;
935: PetscOptionsEnum("-rd_quadrature","Finite element quadrature","",QuadratureTypes,(PetscEnum)rd->quadrature,(PetscEnum*)&rd->quadrature,NULL);
936: }
937: PetscOptionsEnum("-rd_jacobian","Type of finite difference Jacobian","",JacobianTypes,(PetscEnum)rd->jacobian,(PetscEnum*)&rd->jacobian,NULL);
938: switch (rd->initial) {
939: case 1:
940: rd->leftbc = BC_ROBIN;
941: rd->Eapplied = 4 * rd->unit.Joule / PetscPowRealInt(rd->unit.meter,3);
942: rd->L = 1. * rd->unit.meter;
943: rd->beta = 3.0;
944: rd->gamma = 3.0;
945: rd->final_time = 3 * second;
946: break;
947: case 2:
948: rd->leftbc = BC_NEUMANN;
949: rd->Eapplied = 0.;
950: rd->L = 1. * rd->unit.meter;
951: rd->beta = 3.0;
952: rd->gamma = 3.0;
953: rd->final_time = 1 * second;
954: break;
955: case 3:
956: rd->leftbc = BC_ROBIN;
957: rd->Eapplied = 7.503e6 * rd->unit.Joule / PetscPowRealInt(rd->unit.meter,3);
958: rd->L = 5. * rd->unit.meter;
959: rd->beta = 3.5;
960: rd->gamma = 3.5;
961: rd->final_time = 20e-9 * second;
962: break;
963: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Initial %D",rd->initial);
964: }
965: PetscOptionsEnum("-rd_leftbc","Left boundary condition","",BCTypes,(PetscEnum)rd->leftbc,(PetscEnum*)&rd->leftbc,NULL);
966: PetscOptionsReal("-rd_E_applied","Radiation flux at left end of domain","",rd->Eapplied,&rd->Eapplied,NULL);
967: PetscOptionsReal("-rd_beta","Thermal exponent for photon absorption","",rd->beta,&rd->beta,NULL);
968: PetscOptionsReal("-rd_gamma","Thermal exponent for diffusion coefficient","",rd->gamma,&rd->gamma,NULL);
969: PetscOptionsBool("-rd_view_draw","Draw final solution","",rd->view_draw,&rd->view_draw,NULL);
970: PetscOptionsBool("-rd_endpoint","Discretize using endpoints (like trapezoid rule) instead of midpoint","",rd->endpoint,&rd->endpoint,NULL);
971: PetscOptionsBool("-rd_bcmidpoint","Impose the boundary condition at the midpoint (Theta) of the interval","",rd->bcmidpoint,&rd->bcmidpoint,NULL);
972: PetscOptionsBool("-rd_bclimit","Limit diffusion coefficient in definition of Robin boundary condition","",rd->bclimit,&rd->bclimit,NULL);
973: PetscOptionsBool("-rd_test_diff","Test differentiation in constitutive relations","",rd->test_diff,&rd->test_diff,NULL);
974: PetscOptionsString("-rd_view_binary","File name to hold final solution","",rd->view_binary,rd->view_binary,sizeof(rd->view_binary),NULL);
975: }
976: PetscOptionsEnd();
978: switch (rd->initial) {
979: case 1:
980: case 2:
981: rd->rho = 1.;
982: rd->c = 1.;
983: rd->K_R = 1.;
984: rd->K_p = 1.;
985: rd->sigma_b = 0.25;
986: rd->MaterialEnergy = RDMaterialEnergy_Reduced;
987: break;
988: case 3:
989: /* Table 2 */
990: rd->rho = 1.17e-3 * kilogram / (meter*meter*meter); /* density */
991: rd->K_R = 7.44e18 * PetscPowRealInt(meter,5) * PetscPowReal(Kelvin,3.5) * PetscPowRealInt(kilogram,-2); /* */
992: rd->K_p = 2.33e20 * PetscPowRealInt(meter,5) * PetscPowReal(Kelvin,3.5) * PetscPowRealInt(kilogram,-2); /* */
993: rd->I_H = 2.179e-18 * Joule; /* Hydrogen ionization potential */
994: rd->m_p = 1.673e-27 * kilogram; /* proton mass */
995: rd->m_e = 9.109e-31 * kilogram; /* electron mass */
996: rd->h = 6.626e-34 * Joule * second; /* Planck's constant */
997: rd->k = 1.381e-23 * Joule / Kelvin; /* Boltzman constant */
998: rd->c = 3.00e8 * meter / second; /* speed of light */
999: rd->sigma_b = 5.67e-8 * Watt * PetscPowRealInt(meter,-2) * PetscPowRealInt(Kelvin,-4); /* Stefan-Boltzman constant */
1000: rd->MaterialEnergy = RDMaterialEnergy_Saha;
1001: break;
1002: }
1004: DMDACreate1d(comm,DM_BOUNDARY_NONE,20,sizeof(RDNode)/sizeof(PetscScalar),1,NULL,&rd->da);
1005: DMSetFromOptions(rd->da);
1006: DMSetUp(rd->da);
1007: DMDASetFieldName(rd->da,0,"E");
1008: DMDASetFieldName(rd->da,1,"T");
1009: DMDASetUniformCoordinates(rd->da,0.,1.,0.,0.,0.,0.);
1011: *inrd = rd;
1012: return 0;
1013: }
1015: int main(int argc, char *argv[])
1016: {
1017: RD rd;
1018: TS ts;
1019: SNES snes;
1020: Vec X;
1021: Mat A,B;
1022: PetscInt steps;
1023: PetscReal ftime;
1025: PetscInitialize(&argc,&argv,0,help);
1026: RDCreate(PETSC_COMM_WORLD,&rd);
1027: DMCreateGlobalVector(rd->da,&X);
1028: DMSetMatType(rd->da,MATAIJ);
1029: DMCreateMatrix(rd->da,&B);
1030: RDInitialState(rd,X);
1032: TSCreate(PETSC_COMM_WORLD,&ts);
1033: TSSetProblemType(ts,TS_NONLINEAR);
1034: TSSetType(ts,TSTHETA);
1035: TSSetDM(ts,rd->da);
1036: switch (rd->discretization) {
1037: case DISCRETIZATION_FD:
1038: TSSetIFunction(ts,NULL,RDIFunction_FD,rd);
1039: if (rd->jacobian == JACOBIAN_ANALYTIC) TSSetIJacobian(ts,B,B,RDIJacobian_FD,rd);
1040: break;
1041: case DISCRETIZATION_FE:
1042: TSSetIFunction(ts,NULL,RDIFunction_FE,rd);
1043: if (rd->jacobian == JACOBIAN_ANALYTIC) TSSetIJacobian(ts,B,B,RDIJacobian_FE,rd);
1044: break;
1045: }
1046: TSSetMaxTime(ts,rd->final_time);
1047: TSSetTimeStep(ts,1e-3);
1048: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
1049: TSSetFromOptions(ts);
1051: A = B;
1052: TSGetSNES(ts,&snes);
1053: switch (rd->jacobian) {
1054: case JACOBIAN_ANALYTIC:
1055: break;
1056: case JACOBIAN_MATRIXFREE:
1057: break;
1058: case JACOBIAN_FD_COLORING: {
1059: SNESSetJacobian(snes,A,B,SNESComputeJacobianDefaultColor,0);
1060: } break;
1061: case JACOBIAN_FD_FULL:
1062: SNESSetJacobian(snes,A,B,SNESComputeJacobianDefault,ts);
1063: break;
1064: }
1066: if (rd->test_diff) {
1067: RDTestDifferentiation(rd);
1068: }
1069: TSSolve(ts,X);
1070: TSGetSolveTime(ts,&ftime);
1071: TSGetStepNumber(ts,&steps);
1072: PetscPrintf(PETSC_COMM_WORLD,"Steps %D final time %g\n",steps,(double)ftime);
1073: if (rd->view_draw) {
1074: RDView(rd,X,PETSC_VIEWER_DRAW_WORLD);
1075: }
1076: if (rd->view_binary[0]) {
1077: PetscViewer viewer;
1078: PetscViewerBinaryOpen(PETSC_COMM_WORLD,rd->view_binary,FILE_MODE_WRITE,&viewer);
1079: RDView(rd,X,viewer);
1080: PetscViewerDestroy(&viewer);
1081: }
1082: VecDestroy(&X);
1083: MatDestroy(&B);
1084: RDDestroy(&rd);
1085: TSDestroy(&ts);
1086: PetscFinalize();
1087: return 0;
1088: }
1089: /*TEST
1091: test:
1092: 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
1093: requires: !single
1095: test:
1096: suffix: 2
1097: 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
1098: requires: !single
1100: test:
1101: suffix: 3
1102: nsize: 2
1103: 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
1104: requires: !single
1106: TEST*/