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