Actual source code: ex10.c

petsc-3.14.6 2021-03-30
Report Typos and Errors
  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: {

 70:   DMDestroy(&(*rd)->da);
 71:   PetscFree(*rd);
 72:   return(0);
 73: }

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

 85: /* Solves a quadratic equation while propagating tangents */
 86: static void QuadraticSolve(PetscScalar a,PetscScalar a_t,PetscScalar b,PetscScalar b_t,PetscScalar c,PetscScalar c_t,PetscScalar *x,PetscScalar *x_t)
 87: {
 88:   PetscScalar
 89:     disc   = b*b - 4.*a*c,
 90:     disc_t = 2.*b*b_t - 4.*a_t*c - 4.*a*c_t,
 91:     num    = -b + PetscSqrtScalar(disc), /* choose positive sign */
 92:     num_t  = -b_t + 0.5/PetscSqrtScalar(disc)*disc_t,
 93:     den    = 2.*a,
 94:     den_t  = 2.*a_t;
 95:   *x   = num/den;
 96:   *x_t = (num_t*den - num*den_t) / PetscSqr(den);
 97: }

 99: /* The primary model presented in the paper */
100: static void RDMaterialEnergy_Saha(RD rd,const RDNode *n,PetscScalar *inEm,RDNode *dEm)
101: {
102:   PetscScalar Em,alpha,alpha_t,
103:               T     = n->T,
104:               T_t   = 1.,
105:               chi   = rd->I_H / (rd->k * T),
106:               chi_t = -chi / T * T_t,
107:               a     = 1.,
108:               a_t   = 0,
109:               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 */
110:               b_t   = -b*chi_t + 1.5*b/chi*chi_t,
111:               c     = -b,
112:               c_t   = -b_t;
113:   QuadraticSolve(a,a_t,b,b_t,c,c_t,&alpha,&alpha_t);       /* Solve Eq 7 for alpha */
114:   Em = rd->k * T / rd->m_p * (1.5*(1.+alpha) + alpha*chi); /* Eq 6 */
115:   if (inEm) *inEm = Em;
116:   if (dEm) {
117:     dEm->E = 0;
118:     dEm->T = Em / T * T_t + rd->k * T / rd->m_p * (1.5*alpha_t + alpha_t*chi + alpha*chi_t);
119:   }
120: }
121: /* Reduced ionization model, Eq 30 */
122: static void RDMaterialEnergy_Reduced(RD rd,const RDNode *n,PetscScalar *Em,RDNode *dEm)
123: {
124:   PetscScalar alpha,alpha_t,
125:               T     = n->T,
126:               T_t   = 1.,
127:               chi   = -0.3 / T,
128:               chi_t = -chi / T * T_t,
129:               a     = 1.,
130:               a_t   = 0.,
131:               b     = PetscExpScalar(chi),
132:               b_t   = b*chi_t,
133:               c     = -b,
134:               c_t   = -b_t;
135:   QuadraticSolve(a,a_t,b,b_t,c,c_t,&alpha,&alpha_t);
136:   if (Em) *Em = (1.+alpha)*T + 0.3*alpha;
137:   if (dEm) {
138:     dEm->E = 0;
139:     dEm->T = alpha_t*T + (1.+alpha)*T_t + 0.3*alpha_t;
140:   }
141: }

143: /* Eq 5 */
144: static void RDSigma_R(RD rd,RDNode *n,PetscScalar *sigma_R,RDNode *dsigma_R)
145: {
146:   *sigma_R    = rd->K_R * rd->rho * PetscPowScalar(n->T,-rd->gamma);
147:   dsigma_R->E = 0;
148:   dsigma_R->T = -rd->gamma * (*sigma_R) / n->T;
149: }

151: /* Eq 4 */
152: static void RDDiffusionCoefficient(RD rd,PetscBool limit,RDNode *n,RDNode *nx,PetscScalar *D_R,RDNode *dD_R,RDNode *dxD_R)
153: {
154:   PetscScalar sigma_R,denom;
155:   RDNode      dsigma_R,ddenom,dxdenom;

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

174: static PetscErrorCode RDStateView(RD rd,Vec X,Vec Xdot,Vec F)
175: {
177:   DMDALocalInfo  info;
178:   PetscInt       i;
179:   const RDNode   *x,*xdot,*f;
180:   MPI_Comm       comm;

183:   PetscObjectGetComm((PetscObject)rd->da,&comm);
184:   DMDAGetLocalInfo(rd->da,&info);
185:   DMDAVecGetArrayRead(rd->da,X,(void*)&x);
186:   DMDAVecGetArrayRead(rd->da,Xdot,(void*)&xdot);
187:   DMDAVecGetArrayRead(rd->da,F,(void*)&f);
188:   for (i=info.xs; i<info.xs+info.xm; i++) {
189:     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),
190:                                    PetscRealPart(xdot[i].E),PetscRealPart(xdot[i].T), PetscRealPart(f[i].E),PetscRealPart(f[i].T));
191:   }
192:   DMDAVecRestoreArrayRead(rd->da,X,(void*)&x);
193:   DMDAVecRestoreArrayRead(rd->da,Xdot,(void*)&xdot);
194:   DMDAVecRestoreArrayRead(rd->da,F,(void*)&f);
195:   PetscSynchronizedFlush(comm,PETSC_STDOUT);
196:   return(0);
197: }

199: static PetscScalar RDRadiation(RD rd,const RDNode *n,RDNode *dn)
200: {
201:   PetscScalar sigma_p   = rd->K_p * rd->rho * PetscPowScalar(n->T,-rd->beta),
202:               sigma_p_T = -rd->beta * sigma_p / n->T,
203:               tmp       = 4.* rd->sigma_b*PetscSqr(PetscSqr(n->T)) / rd->c - n->E,
204:               tmp_E     = -1.,
205:               tmp_T     = 4. * rd->sigma_b * 4 * n->T*(PetscSqr(n->T)) / rd->c,
206:               rad       = sigma_p * rd->c * rd->rho * tmp,
207:               rad_E     = sigma_p * rd->c * rd->rho * tmp_E,
208:               rad_T     = rd->c * rd->rho * (sigma_p_T * tmp + sigma_p * tmp_T);
209:   if (dn) {
210:     dn->E = rad_E;
211:     dn->T = rad_T;
212:   }
213:   return rad;
214: }

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

222:   n_L.E  = 0.5*(x[i-1].E + x[i].E);
223:   n_L.T  = 0.5*(x[i-1].T + x[i].T);
224:   nx_L.E = (x[i].E - x[i-1].E)/hx;
225:   nx_L.T = (x[i].T - x[i-1].T)/hx;
226:   RDDiffusionCoefficient(rd,PETSC_TRUE,&n_L,&nx_L,&D_L,&dD_L,&dxD_L);
227:   fluxL       = D_L*nx_L.E;
228:   dfluxL[0].E = -ihx*D_L + (0.5*dD_L.E - ihx*dxD_L.E)*nx_L.E;
229:   dfluxL[1].E = +ihx*D_L + (0.5*dD_L.E + ihx*dxD_L.E)*nx_L.E;
230:   dfluxL[0].T = (0.5*dD_L.T - ihx*dxD_L.T)*nx_L.E;
231:   dfluxL[1].T = (0.5*dD_L.T + ihx*dxD_L.T)*nx_L.E;

233:   n_R.E  = 0.5*(x[i].E + x[i+1].E);
234:   n_R.T  = 0.5*(x[i].T + x[i+1].T);
235:   nx_R.E = (x[i+1].E - x[i].E)/hx;
236:   nx_R.T = (x[i+1].T - x[i].T)/hx;
237:   RDDiffusionCoefficient(rd,PETSC_TRUE,&n_R,&nx_R,&D_R,&dD_R,&dxD_R);
238:   fluxR       = D_R*nx_R.E;
239:   dfluxR[0].E = -ihx*D_R + (0.5*dD_R.E - ihx*dxD_R.E)*nx_R.E;
240:   dfluxR[1].E = +ihx*D_R + (0.5*dD_R.E + ihx*dxD_R.E)*nx_R.E;
241:   dfluxR[0].T = (0.5*dD_R.T - ihx*dxD_R.T)*nx_R.E;
242:   dfluxR[1].T = (0.5*dD_R.T + ihx*dxD_R.T)*nx_R.E;

244:   if (d) {
245:     d[0].E = -ihx*dfluxL[0].E;
246:     d[0].T = -ihx*dfluxL[0].T;
247:     d[1].E =  ihx*(dfluxR[0].E - dfluxL[1].E);
248:     d[1].T =  ihx*(dfluxR[0].T - dfluxL[1].T);
249:     d[2].E =  ihx*dfluxR[1].E;
250:     d[2].T =  ihx*dfluxR[1].T;
251:   }
252:   return ihx*(fluxR - fluxL);
253: }

255: 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)
256: {
258:   PetscBool      istheta;

261:   DMGetLocalVector(rd->da,X0loc);
262:   DMGetLocalVector(rd->da,Xloc);
263:   DMGetLocalVector(rd->da,Xloc_t);

265:   DMGlobalToLocalBegin(rd->da,X,INSERT_VALUES,*Xloc);
266:   DMGlobalToLocalEnd(rd->da,X,INSERT_VALUES,*Xloc);
267:   DMGlobalToLocalBegin(rd->da,Xdot,INSERT_VALUES,*Xloc_t);
268:   DMGlobalToLocalEnd(rd->da,Xdot,INSERT_VALUES,*Xloc_t);

270:   /*
271:     The following is a hack to subvert TSTHETA which is like an implicit midpoint method to behave more like a trapezoid
272:     rule.  These methods have equivalent linear stability, but the nonlinear stability is somewhat different.  The
273:     radiation system is inconvenient to write in explicit form because the ionization model is "on the left".
274:    */
275:   PetscObjectTypeCompare((PetscObject)ts,TSTHETA,&istheta);
276:   if (istheta && rd->endpoint) {
277:     TSThetaGetTheta(ts,Theta);
278:   } else *Theta = 1.;

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

286:   DMDAVecGetArray(rd->da,*X0loc,x0);
287:   DMDAVecGetArray(rd->da,*Xloc,x);
288:   DMDAVecGetArray(rd->da,*Xloc_t,xdot);
289:   return(0);
290: }

292: static PetscErrorCode RDRestoreLocalArrays(RD rd,Vec *X0loc,RDNode **x0,Vec *Xloc,RDNode **x,Vec *Xloc_t,RDNode **xdot)
293: {

297:   DMDAVecRestoreArray(rd->da,*X0loc,x0);
298:   DMDAVecRestoreArray(rd->da,*Xloc,x);
299:   DMDAVecRestoreArray(rd->da,*Xloc_t,xdot);
300:   DMRestoreLocalVector(rd->da,X0loc);
301:   DMRestoreLocalVector(rd->da,Xloc);
302:   DMRestoreLocalVector(rd->da,Xloc_t);
303:   return(0);
304: }

306: static PetscErrorCode RDCheckDomain_Private(RD rd,TS ts,Vec X,PetscBool  *in)
307: {
309:   PetscInt       minloc;
310:   PetscReal      min;

313:   VecMin(X,&minloc,&min);
314:   if (min < 0) {
315:     SNES snes;
316:     *in  = PETSC_FALSE;
317:     TSGetSNES(ts,&snes);
318:     SNESSetFunctionDomainError(snes);
319:     PetscInfo3(ts,"Domain violation at %D field %D value %g\n",minloc/2,minloc%2,(double)min);
320:   } else *in = PETSC_TRUE;
321:   return(0);
322: }

324: /* Energy and temperature must remain positive */
325: #define RDCheckDomain(rd,ts,X) do {                                \
326:     PetscErrorCode _ierr;                                          \
327:     PetscBool      _in;                                                \
328:     _RDCheckDomain_Private(rd,ts,X,&_in);CHKERRQ(_ierr);    \
329:     if (!_in) return(0);                              \
330:   } while (0)

332: static PetscErrorCode RDIFunction_FD(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
333: {
335:   RD             rd = (RD)ctx;
336:   RDNode         *x,*x0,*xdot,*f;
337:   Vec            X0loc,Xloc,Xloc_t;
338:   PetscReal      hx,Theta,dt;
339:   DMDALocalInfo  info;
340:   PetscInt       i;

343:   RDCheckDomain(rd,ts,X);
344:   RDGetLocalArrays(rd,ts,X,Xdot,&Theta,&dt,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
345:   DMDAVecGetArray(rd->da,F,&f);
346:   DMDAGetLocalInfo(rd->da,&info);
347:   VecZeroEntries(F);

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

351:   for (i=info.xs; i<info.xs+info.xm; i++) {
352:     PetscReal   rho = rd->rho;
353:     PetscScalar Em_t,rad;

355:     rad = (1.-Theta)*RDRadiation(rd,&x0[i],0) + Theta*RDRadiation(rd,&x[i],0);
356:     if (rd->endpoint) {
357:       PetscScalar Em0,Em1;
358:       RDMaterialEnergy(rd,&x0[i],&Em0,NULL);
359:       RDMaterialEnergy(rd,&x[i],&Em1,NULL);
360:       Em_t = (Em1 - Em0) / dt;
361:     } else {
362:       RDNode dEm;
363:       RDMaterialEnergy(rd,&x[i],NULL,&dEm);
364:       Em_t = dEm.E * xdot[i].E + dEm.T * xdot[i].T;
365:     }
366:     /* Residuals are multiplied by the volume element (hx).  */
367:     /* The temperature equation does not have boundary conditions */
368:     f[i].T = hx*(rho*Em_t + rad);

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

374:       n.E  =  (1.-bcTheta)*x0[0].E + bcTheta*x[0].E;
375:       n.T  =  (1.-bcTheta)*x0[0].T + bcTheta*x[0].T;
376:       nx.E = ((1.-bcTheta)*(x0[1].E-x0[0].E) + bcTheta*(x[1].E-x[0].E))/hx;
377:       nx.T = ((1.-bcTheta)*(x0[1].T-x0[0].T) + bcTheta*(x[1].T-x[0].T))/hx;
378:       switch (rd->leftbc) {
379:       case BC_ROBIN:
380:         RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D_R,0,0);
381:         f[0].E = hx*(n.E - 2. * D_R * nx.E - rd->Eapplied);
382:         break;
383:       case BC_NEUMANN:
384:         f[0].E = x[1].E - x[0].E;
385:         break;
386:       default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D",rd->initial);
387:       }
388:     } else if (i == info.mx-1) { /* Right boundary */
389:       f[i].E = x[i].E - x[i-1].E; /* Homogeneous Neumann */
390:     } else {
391:       PetscScalar diff = (1.-Theta)*RDDiffusion(rd,hx,x0,i,0) + Theta*RDDiffusion(rd,hx,x,i,0);
392:       f[i].E = hx*(xdot[i].E - diff - rad);
393:     }
394:   }
395:   RDRestoreLocalArrays(rd,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
396:   DMDAVecRestoreArray(rd->da,F,&f);
397:   if (rd->monitor_residual) {RDStateView(rd,X,Xdot,F);}
398:   return(0);
399: }

401: static PetscErrorCode RDIJacobian_FD(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
402: {
404:   RD             rd = (RD)ctx;
405:   RDNode         *x,*x0,*xdot;
406:   Vec            X0loc,Xloc,Xloc_t;
407:   PetscReal      hx,Theta,dt;
408:   DMDALocalInfo  info;
409:   PetscInt       i;

412:   RDCheckDomain(rd,ts,X);
413:   RDGetLocalArrays(rd,ts,X,Xdot,&Theta,&dt,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
414:   DMDAGetLocalInfo(rd->da,&info);
415:   hx   = rd->L / (info.mx-1);
416:   MatZeroEntries(B);

418:   for (i=info.xs; i<info.xs+info.xm; i++) {
419:     PetscInt                  col[3];
420:     PetscReal                 rho = rd->rho;
421:     PetscScalar /*Em_t,rad,*/ K[2][6];
422:     RDNode                    dEm_t,drad;

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

426:     if (rd->endpoint) {
427:       PetscScalar Em0,Em1;
428:       RDNode      dEm1;
429:       RDMaterialEnergy(rd,&x0[i],&Em0,NULL);
430:       RDMaterialEnergy(rd,&x[i],&Em1,&dEm1);
431:       /*Em_t = (Em1 - Em0) / (Theta*dt);*/
432:       dEm_t.E = dEm1.E / (Theta*dt);
433:       dEm_t.T = dEm1.T / (Theta*dt);
434:     } else {
435:       const PetscScalar epsilon = x[i].T * PETSC_SQRT_MACHINE_EPSILON;
436:       RDNode            n1;
437:       RDNode            dEm,dEm1;
438:       PetscScalar       Em_TT;

440:       n1.E = x[i].E;
441:       n1.T = x[i].T+epsilon;
442:       RDMaterialEnergy(rd,&x[i],NULL,&dEm);
443:       RDMaterialEnergy(rd,&n1,NULL,&dEm1);
444:       /* The Jacobian needs another derivative.  We finite difference here instead of
445:        * propagating second derivatives through the ionization model. */
446:       Em_TT = (dEm1.T - dEm.T) / epsilon;
447:       /*Em_t = dEm.E * xdot[i].E + dEm.T * xdot[i].T;*/
448:       dEm_t.E = dEm.E * a;
449:       dEm_t.T = dEm.T * a + Em_TT * xdot[i].T;
450:     }

452:     PetscMemzero(K,sizeof(K));
453:     /* Residuals are multiplied by the volume element (hx).  */
454:     if (i == 0) {
455:       PetscScalar D,bcTheta = rd->bcmidpoint ? Theta : 1.;
456:       RDNode      n, nx;
457:       RDNode      dD,dxD;

459:       n.E  = (1.-bcTheta)*x0[0].E + bcTheta*x[0].E;
460:       n.T  = (1.-bcTheta)*x0[0].T + bcTheta*x[0].T;
461:       nx.E = ((1.-bcTheta)*(x0[1].E-x0[0].E) + bcTheta*(x[1].E-x[0].E))/hx;
462:       nx.T = ((1.-bcTheta)*(x0[1].T-x0[0].T) + bcTheta*(x[1].T-x[0].T))/hx;
463:       switch (rd->leftbc) {
464:       case BC_ROBIN:
465:         RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D,&dD,&dxD);
466:         K[0][1*2+0] = (bcTheta/Theta)*hx*(1. -2.*D*(-1./hx) - 2.*nx.E*dD.E + 2.*nx.E*dxD.E/hx);
467:         K[0][1*2+1] = (bcTheta/Theta)*hx*(-2.*nx.E*dD.T);
468:         K[0][2*2+0] = (bcTheta/Theta)*hx*(-2.*D*(1./hx) - 2.*nx.E*dD.E - 2.*nx.E*dxD.E/hx);
469:         break;
470:       case BC_NEUMANN:
471:         K[0][1*2+0] = -1./Theta;
472:         K[0][2*2+0] = 1./Theta;
473:         break;
474:       default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D",rd->initial);
475:       }
476:     } else if (i == info.mx-1) {
477:       K[0][0*2+0] = -1./Theta;
478:       K[0][1*2+0] = 1./Theta;
479:     } else {
480:       /*PetscScalar diff;*/
481:       RDNode ddiff[3];
482:       /*diff = (1.-Theta)*RDDiffusion(rd,hx,x0,i,0) + Theta* */ RDDiffusion(rd,hx,x,i,ddiff);
483:       K[0][0*2+0] = -hx*ddiff[0].E;
484:       K[0][0*2+1] = -hx*ddiff[0].T;
485:       K[0][1*2+0] = hx*(a - ddiff[1].E - drad.E);
486:       K[0][1*2+1] = hx*(-ddiff[1].T - drad.T);
487:       K[0][2*2+0] = -hx*ddiff[2].E;
488:       K[0][2*2+1] = -hx*ddiff[2].T;
489:     }

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

494:     col[0] = i-1;
495:     col[1] = i;
496:     col[2] = i+1<info.mx ? i+1 : -1;
497:     MatSetValuesBlocked(B,1,&i,3,col,&K[0][0],INSERT_VALUES);
498:   }
499:   RDRestoreLocalArrays(rd,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
500:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
501:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
502:   if (A != B) {
503:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
504:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
505:   }
506:   return(0);
507: }

509: /* Evaluate interpolants and derivatives at a select quadrature point */
510: static void RDEvaluate(PetscReal interp[][2],PetscReal deriv[][2],PetscInt q,const RDNode x[],PetscInt i,RDNode *n,RDNode *nx)
511: {
512:   PetscInt j;
513:   n->E = 0; n->T = 0; nx->E = 0; nx->T = 0;
514:   for (j=0; j<2; j++) {
515:     n->E  += interp[q][j] * x[i+j].E;
516:     n->T  += interp[q][j] * x[i+j].T;
517:     nx->E += deriv[q][j] * x[i+j].E;
518:     nx->T += deriv[q][j] * x[i+j].T;
519:   }
520: }

522: /*
523:  Various quadrature rules.  The nonlinear terms are non-polynomial so no standard quadrature will be exact.
524: */
525: static PetscErrorCode RDGetQuadrature(RD rd,PetscReal hx,PetscInt *nq,PetscReal weight[],PetscReal interp[][2],PetscReal deriv[][2])
526: {
527:   PetscInt        q,j;
528:   const PetscReal *refweight,(*refinterp)[2],(*refderiv)[2];

531:   switch (rd->quadrature) {
532:   case QUADRATURE_GAUSS1: {
533:     static const PetscReal ww[1] = {1.},ii[1][2] = {{0.5,0.5}},dd[1][2] = {{-1.,1.}};
534:     *nq = 1; refweight = ww; refinterp = ii; refderiv = dd;
535:   } break;
536:   case QUADRATURE_GAUSS2: {
537:     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};
538:     *nq = 2; refweight = ww; refinterp = ii; refderiv = dd;
539:   } break;
540:   case QUADRATURE_GAUSS3: {
541:     static const PetscReal ii[3][2] = {{0.8872983346207417,0.1127016653792583},{0.5,0.5},{0.1127016653792583,0.8872983346207417}},
542:                            dd[3][2] = {{-1,1},{-1,1},{-1,1}},ww[3] = {5./18,8./18,5./18};
543:     *nq = 3; refweight = ww; refinterp = ii; refderiv = dd;
544:   } break;
545:   case QUADRATURE_GAUSS4: {
546:     static const PetscReal ii[][2] = {{0.93056815579702623,0.069431844202973658},
547:                                       {0.66999052179242813,0.33000947820757187},
548:                                       {0.33000947820757187,0.66999052179242813},
549:                                       {0.069431844202973658,0.93056815579702623}},
550:                            dd[][2] = {{-1,1},{-1,1},{-1,1},{-1,1}},ww[] = {0.17392742256872692,0.3260725774312731,0.3260725774312731,0.17392742256872692};

552:     *nq = 4; refweight = ww; refinterp = ii; refderiv = dd;
553:   } break;
554:   case QUADRATURE_LOBATTO2: {
555:     static const PetscReal ii[2][2] = {{1.,0.},{0.,1.}},dd[2][2] = {{-1.,1.},{-1.,1.}},ww[2] = {0.5,0.5};
556:     *nq = 2; refweight = ww; refinterp = ii; refderiv = dd;
557:   } break;
558:   case QUADRATURE_LOBATTO3: {
559:     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};
560:     *nq = 3; refweight = ww; refinterp = ii; refderiv = dd;
561:   } break;
562:   default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown quadrature %d",(int)rd->quadrature);
563:   }

565:   for (q=0; q<*nq; q++) {
566:     weight[q] = refweight[q] * hx;
567:     for (j=0; j<2; j++) {
568:       interp[q][j] = refinterp[q][j];
569:       deriv[q][j]  = refderiv[q][j] / hx;
570:     }
571:   }
572:   return(0);
573: }

575: /*
576:  Finite element version
577: */
578: static PetscErrorCode RDIFunction_FE(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
579: {
581:   RD             rd = (RD)ctx;
582:   RDNode         *x,*x0,*xdot,*f;
583:   Vec            X0loc,Xloc,Xloc_t,Floc;
584:   PetscReal      hx,Theta,dt,weight[5],interp[5][2],deriv[5][2];
585:   DMDALocalInfo  info;
586:   PetscInt       i,j,q,nq;

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

592:   DMGetLocalVector(rd->da,&Floc);
593:   VecZeroEntries(Floc);
594:   DMDAVecGetArray(rd->da,Floc,&f);
595:   DMDAGetLocalInfo(rd->da,&info);

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

601:   for (i=info.xs; i<PetscMin(info.xs+info.xm,info.mx-1); i++) {
602:     for (q=0; q<nq; q++) {
603:       PetscReal   rho = rd->rho;
604:       PetscScalar Em_t,rad,D_R,D0_R;
605:       RDNode      n,n0,nx,n0x,nt,ntx;
606:       RDEvaluate(interp,deriv,q,x,i,&n,&nx);
607:       RDEvaluate(interp,deriv,q,x0,i,&n0,&n0x);
608:       RDEvaluate(interp,deriv,q,xdot,i,&nt,&ntx);

610:       rad = (1.-Theta)*RDRadiation(rd,&n0,0) + Theta*RDRadiation(rd,&n,0);
611:       if (rd->endpoint) {
612:         PetscScalar Em0,Em1;
613:         RDMaterialEnergy(rd,&n0,&Em0,NULL);
614:         RDMaterialEnergy(rd,&n,&Em1,NULL);
615:         Em_t = (Em1 - Em0) / dt;
616:       } else {
617:         RDNode dEm;
618:         RDMaterialEnergy(rd,&n,NULL,&dEm);
619:         Em_t = dEm.E * nt.E + dEm.T * nt.T;
620:       }
621:       RDDiffusionCoefficient(rd,PETSC_TRUE,&n0,&n0x,&D0_R,0,0);
622:       RDDiffusionCoefficient(rd,PETSC_TRUE,&n,&nx,&D_R,0,0);
623:       for (j=0; j<2; j++) {
624:         f[i+j].E += (deriv[q][j] * weight[q] * ((1.-Theta)*D0_R*n0x.E + Theta*D_R*nx.E)
625:                      + interp[q][j] * weight[q] * (nt.E - rad));
626:         f[i+j].T += interp[q][j] * weight[q] * (rho * Em_t + rad);
627:       }
628:     }
629:   }
630:   if (info.xs == 0) {
631:     switch (rd->leftbc) {
632:     case BC_ROBIN: {
633:       PetscScalar D_R,D_R_bc;
634:       PetscReal   ratio,bcTheta = rd->bcmidpoint ? Theta : 1.;
635:       RDNode      n, nx;

637:       n.E  = (1-bcTheta)*x0[0].E + bcTheta*x[0].E;
638:       n.T  = (1-bcTheta)*x0[0].T + bcTheta*x[0].T;
639:       nx.E = (x[1].E-x[0].E)/hx;
640:       nx.T = (x[1].T-x[0].T)/hx;
641:       RDDiffusionCoefficient(rd,PETSC_TRUE,&n,&nx,&D_R,0,0);
642:       RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D_R_bc,0,0);
643:       ratio = PetscRealPart(D_R/D_R_bc);
644:       if (ratio > 1.) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Limited diffusivity is greater than unlimited");
645:       if (ratio < 1e-3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Heavily limited diffusivity");
646:       f[0].E += -ratio*0.5*(rd->Eapplied - n.E);
647:     } break;
648:     case BC_NEUMANN:
649:       /* homogeneous Neumann is the natural condition */
650:       break;
651:     default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D",rd->initial);
652:     }
653:   }

655:   RDRestoreLocalArrays(rd,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
656:   DMDAVecRestoreArray(rd->da,Floc,&f);
657:   VecZeroEntries(F);
658:   DMLocalToGlobalBegin(rd->da,Floc,ADD_VALUES,F);
659:   DMLocalToGlobalEnd(rd->da,Floc,ADD_VALUES,F);
660:   DMRestoreLocalVector(rd->da,&Floc);

662:   if (rd->monitor_residual) {RDStateView(rd,X,Xdot,F);}
663:   return(0);
664: }

666: static PetscErrorCode RDIJacobian_FE(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
667: {
669:   RD             rd = (RD)ctx;
670:   RDNode         *x,*x0,*xdot;
671:   Vec            X0loc,Xloc,Xloc_t;
672:   PetscReal      hx,Theta,dt,weight[5],interp[5][2],deriv[5][2];
673:   DMDALocalInfo  info;
674:   PetscInt       i,j,k,q,nq;
675:   PetscScalar    K[4][4];

678:   RDCheckDomain(rd,ts,X);
679:   RDGetLocalArrays(rd,ts,X,Xdot,&Theta,&dt,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
680:   DMDAGetLocalInfo(rd->da,&info);
681:   hx   = rd->L / (info.mx-1);
682:   RDGetQuadrature(rd,hx,&nq,weight,interp,deriv);
683:   MatZeroEntries(B);
684:   for (i=info.xs; i<PetscMin(info.xs+info.xm,info.mx-1); i++) {
685:     PetscInt rc[2];

687:     rc[0] = i; rc[1] = i+1;
688:     PetscMemzero(K,sizeof(K));
689:     for (q=0; q<nq; q++) {
690:       PetscScalar              D_R;
691:       PETSC_UNUSED PetscScalar rad;
692:       RDNode                   n,nx,nt,ntx,drad,dD_R,dxD_R,dEm;
693:       RDEvaluate(interp,deriv,q,x,i,&n,&nx);
694:       RDEvaluate(interp,deriv,q,xdot,i,&nt,&ntx);
695:       rad = RDRadiation(rd,&n,&drad);
696:       RDDiffusionCoefficient(rd,PETSC_TRUE,&n,&nx,&D_R,&dD_R,&dxD_R);
697:       RDMaterialEnergy(rd,&n,NULL,&dEm);
698:       for (j=0; j<2; j++) {
699:         for (k=0; k<2; k++) {
700:           K[j*2+0][k*2+0] += (+interp[q][j] * weight[q] * (a - drad.E) * interp[q][k]
701:                               + deriv[q][j] * weight[q] * ((D_R + dxD_R.E * nx.E) * deriv[q][k] + dD_R.E * nx.E * interp[q][k]));
702:           K[j*2+0][k*2+1] += (+interp[q][j] * weight[q] * (-drad.T * interp[q][k])
703:                               + deriv[q][j] * weight[q] * (dxD_R.T * deriv[q][k] + dD_R.T * interp[q][k]) * nx.E);
704:           K[j*2+1][k*2+0] +=   interp[q][j] * weight[q] * drad.E * interp[q][k];
705:           K[j*2+1][k*2+1] +=   interp[q][j] * weight[q] * (a * rd->rho * dEm.T + drad.T) * interp[q][k];
706:         }
707:       }
708:     }
709:     MatSetValuesBlocked(B,2,rc,2,rc,&K[0][0],ADD_VALUES);
710:   }
711:   if (info.xs == 0) {
712:     switch (rd->leftbc) {
713:     case BC_ROBIN: {
714:       PetscScalar D_R,D_R_bc;
715:       PetscReal   ratio;
716:       RDNode      n, nx;

718:       n.E  = (1-Theta)*x0[0].E + Theta*x[0].E;
719:       n.T  = (1-Theta)*x0[0].T + Theta*x[0].T;
720:       nx.E = (x[1].E-x[0].E)/hx;
721:       nx.T = (x[1].T-x[0].T)/hx;
722:       RDDiffusionCoefficient(rd,PETSC_TRUE,&n,&nx,&D_R,0,0);
723:       RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D_R_bc,0,0);
724:       ratio = PetscRealPart(D_R/D_R_bc);
725:       MatSetValue(B,0,0,ratio*0.5,ADD_VALUES);
726:     } break;
727:     case BC_NEUMANN:
728:       /* homogeneous Neumann is the natural condition */
729:       break;
730:     default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D",rd->initial);
731:     }
732:   }

734:   RDRestoreLocalArrays(rd,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
735:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
736:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
737:   if (A != B) {
738:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
739:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
740:   }
741:   return(0);
742: }

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

747: static PetscErrorCode RDInitialState(RD rd,Vec X)
748: {
749:   DMDALocalInfo  info;
750:   PetscInt       i;
751:   RDNode         *x;

755:   DMDAGetLocalInfo(rd->da,&info);
756:   DMDAVecGetArray(rd->da,X,&x);
757:   for (i=info.xs; i<info.xs+info.xm; i++) {
758:     PetscReal coord = i*rd->L/(info.mx-1);
759:     switch (rd->initial) {
760:     case 1:
761:       x[i].E = 0.001;
762:       x[i].T = RDRadiationTemperature(rd,x[i].E);
763:       break;
764:     case 2:
765:       x[i].E = 0.001 + 100.*PetscExpReal(-PetscSqr(coord/0.1));
766:       x[i].T = RDRadiationTemperature(rd,x[i].E);
767:       break;
768:     case 3:
769:       x[i].E = 7.56e-2 * rd->unit.Joule / PetscPowScalarInt(rd->unit.meter,3);
770:       x[i].T = RDRadiationTemperature(rd,x[i].E);
771:       break;
772:     default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No initial state %D",rd->initial);
773:     }
774:   }
775:   DMDAVecRestoreArray(rd->da,X,&x);
776:   return(0);
777: }

779: static PetscErrorCode RDView(RD rd,Vec X,PetscViewer viewer)
780: {
782:   Vec            Y;
783:   const RDNode   *x;
784:   PetscScalar    *y;
785:   PetscInt       i,m,M;
786:   const PetscInt *lx;
787:   DM             da;
788:   MPI_Comm       comm;

791:   /*
792:     Create a DMDA (one dof per node, zero stencil width, same layout) to hold Trad
793:     (radiation temperature).  It is not necessary to create a DMDA for this, but this way
794:     output and visualization will have meaningful variable names and correct scales.
795:   */
796:   DMDAGetInfo(rd->da,0, &M,0,0, 0,0,0, 0,0,0,0,0,0);
797:   DMDAGetOwnershipRanges(rd->da,&lx,0,0);
798:   PetscObjectGetComm((PetscObject)rd->da,&comm);
799:   DMDACreate1d(comm,DM_BOUNDARY_NONE,M,1,0,lx,&da);
800:   DMSetFromOptions(da);
801:   DMSetUp(da);
802:   DMDASetUniformCoordinates(da,0.,rd->L,0.,0.,0.,0.);
803:   DMDASetFieldName(da,0,"T_rad");
804:   DMCreateGlobalVector(da,&Y);

806:   /* Compute the radiation temperature from the solution at each node */
807:   VecGetLocalSize(Y,&m);
808:   VecGetArrayRead(X,(const PetscScalar **)&x);
809:   VecGetArray(Y,&y);
810:   for (i=0; i<m; i++) y[i] = RDRadiationTemperature(rd,x[i].E);
811:   VecRestoreArrayRead(X,(const PetscScalar**)&x);
812:   VecRestoreArray(Y,&y);

814:   VecView(Y,viewer);
815:   VecDestroy(&Y);
816:   DMDestroy(&da);
817:   return(0);
818: }

820: static PetscErrorCode RDTestDifferentiation(RD rd)
821: {
822:   MPI_Comm       comm;
824:   RDNode         n,nx;
825:   PetscScalar    epsilon;

828:   PetscObjectGetComm((PetscObject)rd->da,&comm);
829:   epsilon = 1e-8;
830:   {
831:     RDNode      dEm,fdEm;
832:     PetscScalar T0 = 1000.,T1 = T0*(1.+epsilon),Em0,Em1;
833:     n.E = 1.;
834:     n.T = T0;
835:     rd->MaterialEnergy(rd,&n,&Em0,&dEm);
836:     n.E = 1.+epsilon;
837:     n.T = T0;
838:     rd->MaterialEnergy(rd,&n,&Em1,0);
839:     fdEm.E = (Em1-Em0)/epsilon;
840:     n.E = 1.;
841:     n.T = T1;
842:     rd->MaterialEnergy(rd,&n,&Em1,0);
843:     fdEm.T = (Em1-Em0)/(T0*epsilon);
844:     PetscPrintf(comm,"dEm {%g,%g}, fdEm {%g,%g}, diff {%g,%g}\n",(double)PetscRealPart(dEm.E),(double)PetscRealPart(dEm.T),
845:                        (double)PetscRealPart(fdEm.E),(double)PetscRealPart(fdEm.T),(double)PetscRealPart(dEm.E-fdEm.E),(double)PetscRealPart(dEm.T-fdEm.T));
846:   }
847:   {
848:     PetscScalar D0,D;
849:     RDNode      dD,dxD,fdD,fdxD;
850:     n.E = 1.;          n.T = 1.;           nx.E = 1.;          n.T = 1.;
851:     RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D0,&dD,&dxD);
852:     n.E = 1.+epsilon;  n.T = 1.;           nx.E = 1.;          n.T = 1.;
853:     RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D,0,0); fdD.E = (D-D0)/epsilon;
854:     n.E = 1;           n.T = 1.+epsilon;   nx.E = 1.;          n.T = 1.;
855:     RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D,0,0); fdD.T = (D-D0)/epsilon;
856:     n.E = 1;           n.T = 1.;           nx.E = 1.+epsilon;  n.T = 1.;
857:     RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D,0,0); fdxD.E = (D-D0)/epsilon;
858:     n.E = 1;           n.T = 1.;           nx.E = 1.;          n.T = 1.+epsilon;
859:     RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D,0,0); fdxD.T = (D-D0)/epsilon;
860:     PetscPrintf(comm,"dD {%g,%g}, fdD {%g,%g}, diff {%g,%g}\n",(double)PetscRealPart(dD.E),(double)PetscRealPart(dD.T),
861:                        (double)PetscRealPart(fdD.E),(double)PetscRealPart(fdD.T),(double)PetscRealPart(dD.E-fdD.E),(double)PetscRealPart(dD.T-fdD.T));
862:     PetscPrintf(comm,"dxD {%g,%g}, fdxD {%g,%g}, diffx {%g,%g}\n",(double)PetscRealPart(dxD.E),(double)PetscRealPart(dxD.T),
863:                        (double)PetscRealPart(fdxD.E),(double)PetscRealPart(fdxD.T),(double)PetscRealPart(dxD.E-fdxD.E),(double)PetscRealPart(dxD.T-fdxD.T));
864:   }
865:   {
866:     PetscInt    i;
867:     PetscReal   hx = 1.;
868:     PetscScalar a0;
869:     RDNode      n0[3],n1[3],d[3],fd[3];

871:     n0[0].E = 1.;
872:     n0[0].T = 1.;
873:     n0[1].E = 5.;
874:     n0[1].T = 3.;
875:     n0[2].E = 4.;
876:     n0[2].T = 2.;
877:     a0 = RDDiffusion(rd,hx,n0,1,d);
878:     for (i=0; i<3; i++) {
879:       PetscMemcpy(n1,n0,sizeof(n0)); n1[i].E += epsilon;
880:       fd[i].E = (RDDiffusion(rd,hx,n1,1,0)-a0)/epsilon;
881:       PetscMemcpy(n1,n0,sizeof(n0)); n1[i].T += epsilon;
882:       fd[i].T = (RDDiffusion(rd,hx,n1,1,0)-a0)/epsilon;
883:       PetscPrintf(comm,"ddiff[%D] {%g,%g}, fd {%g %g}, diff {%g,%g}\n",i,(double)PetscRealPart(d[i].E),(double)PetscRealPart(d[i].T),
884:                             (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));
885:     }
886:   }
887:   {
888:     PetscScalar rad0,rad;
889:     RDNode      drad,fdrad;
890:     n.E  = 1.;         n.T = 1.;
891:     rad0 = RDRadiation(rd,&n,&drad);
892:     n.E  = 1.+epsilon; n.T = 1.;
893:     rad  = RDRadiation(rd,&n,0); fdrad.E = (rad-rad0)/epsilon;
894:     n.E  = 1.;         n.T = 1.+epsilon;
895:     rad  = RDRadiation(rd,&n,0); fdrad.T = (rad-rad0)/epsilon;
896:     PetscPrintf(comm,"drad {%g,%g}, fdrad {%g,%g}, diff {%g,%g}\n",(double)PetscRealPart(drad.E),(double)PetscRealPart(drad.T),
897:                        (double)PetscRealPart(fdrad.E),(double)PetscRealPart(fdrad.T),(double)PetscRealPart(drad.E-drad.E),(double)PetscRealPart(drad.T-fdrad.T));
898:   }
899:   return(0);
900: }

902: static PetscErrorCode RDCreate(MPI_Comm comm,RD *inrd)
903: {
905:   RD             rd;
906:   PetscReal      meter=0,kilogram=0,second=0,Kelvin=0,Joule=0,Watt=0;

909:   *inrd = 0;
910:   PetscNew(&rd);

912:   PetscOptionsBegin(comm,NULL,"Options for nonequilibrium radiation-diffusion with RD ionization",NULL);
913:   {
914:     rd->initial = 1;
915:     PetscOptionsInt("-rd_initial","Initial condition (1=Marshak, 2=Blast, 3=Marshak+)","",rd->initial,&rd->initial,0);
916:     switch (rd->initial) {
917:     case 1:
918:     case 2:
919:       rd->unit.kilogram = 1.;
920:       rd->unit.meter    = 1.;
921:       rd->unit.second   = 1.;
922:       rd->unit.Kelvin   = 1.;
923:       break;
924:     case 3:
925:       rd->unit.kilogram = 1.e12;
926:       rd->unit.meter    = 1.;
927:       rd->unit.second   = 1.e9;
928:       rd->unit.Kelvin   = 1.;
929:       break;
930:     default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown initial condition %d",rd->initial);
931:     }
932:     /* Fundamental units */
933:     PetscOptionsReal("-rd_unit_meter","Length of 1 meter in nondimensional units","",rd->unit.meter,&rd->unit.meter,0);
934:     PetscOptionsReal("-rd_unit_kilogram","Mass of 1 kilogram in nondimensional units","",rd->unit.kilogram,&rd->unit.kilogram,0);
935:     PetscOptionsReal("-rd_unit_second","Time of a second in nondimensional units","",rd->unit.second,&rd->unit.second,0);
936:     PetscOptionsReal("-rd_unit_Kelvin","Temperature of a Kelvin in nondimensional units","",rd->unit.Kelvin,&rd->unit.Kelvin,0);
937:     /* Derived units */
938:     rd->unit.Joule = rd->unit.kilogram*PetscSqr(rd->unit.meter/rd->unit.second);
939:     rd->unit.Watt  = rd->unit.Joule/rd->unit.second;
940:     /* Local aliases */
941:     meter    = rd->unit.meter;
942:     kilogram = rd->unit.kilogram;
943:     second   = rd->unit.second;
944:     Kelvin   = rd->unit.Kelvin;
945:     Joule    = rd->unit.Joule;
946:     Watt     = rd->unit.Watt;

948:     PetscOptionsBool("-rd_monitor_residual","Display residuals every time they are evaluated","",rd->monitor_residual,&rd->monitor_residual,NULL);
949:     PetscOptionsEnum("-rd_discretization","Discretization type","",DiscretizationTypes,(PetscEnum)rd->discretization,(PetscEnum*)&rd->discretization,NULL);
950:     if (rd->discretization == DISCRETIZATION_FE) {
951:       rd->quadrature = QUADRATURE_GAUSS2;
952:       PetscOptionsEnum("-rd_quadrature","Finite element quadrature","",QuadratureTypes,(PetscEnum)rd->quadrature,(PetscEnum*)&rd->quadrature,NULL);
953:     }
954:     PetscOptionsEnum("-rd_jacobian","Type of finite difference Jacobian","",JacobianTypes,(PetscEnum)rd->jacobian,(PetscEnum*)&rd->jacobian,NULL);
955:     switch (rd->initial) {
956:     case 1:
957:       rd->leftbc     = BC_ROBIN;
958:       rd->Eapplied   = 4 * rd->unit.Joule / PetscPowRealInt(rd->unit.meter,3);
959:       rd->L          = 1. * rd->unit.meter;
960:       rd->beta       = 3.0;
961:       rd->gamma      = 3.0;
962:       rd->final_time = 3 * second;
963:       break;
964:     case 2:
965:       rd->leftbc     = BC_NEUMANN;
966:       rd->Eapplied   = 0.;
967:       rd->L          = 1. * rd->unit.meter;
968:       rd->beta       = 3.0;
969:       rd->gamma      = 3.0;
970:       rd->final_time = 1 * second;
971:       break;
972:     case 3:
973:       rd->leftbc     = BC_ROBIN;
974:       rd->Eapplied   = 7.503e6 * rd->unit.Joule / PetscPowRealInt(rd->unit.meter,3);
975:       rd->L          = 5. * rd->unit.meter;
976:       rd->beta       = 3.5;
977:       rd->gamma      = 3.5;
978:       rd->final_time = 20e-9 * second;
979:       break;
980:     default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Initial %D",rd->initial);
981:     }
982:     PetscOptionsEnum("-rd_leftbc","Left boundary condition","",BCTypes,(PetscEnum)rd->leftbc,(PetscEnum*)&rd->leftbc,NULL);
983:     PetscOptionsReal("-rd_E_applied","Radiation flux at left end of domain","",rd->Eapplied,&rd->Eapplied,NULL);
984:     PetscOptionsReal("-rd_beta","Thermal exponent for photon absorption","",rd->beta,&rd->beta,NULL);
985:     PetscOptionsReal("-rd_gamma","Thermal exponent for diffusion coefficient","",rd->gamma,&rd->gamma,NULL);
986:     PetscOptionsBool("-rd_view_draw","Draw final solution","",rd->view_draw,&rd->view_draw,NULL);
987:     PetscOptionsBool("-rd_endpoint","Discretize using endpoints (like trapezoid rule) instead of midpoint","",rd->endpoint,&rd->endpoint,NULL);
988:     PetscOptionsBool("-rd_bcmidpoint","Impose the boundary condition at the midpoint (Theta) of the interval","",rd->bcmidpoint,&rd->bcmidpoint,NULL);
989:     PetscOptionsBool("-rd_bclimit","Limit diffusion coefficient in definition of Robin boundary condition","",rd->bclimit,&rd->bclimit,NULL);
990:     PetscOptionsBool("-rd_test_diff","Test differentiation in constitutive relations","",rd->test_diff,&rd->test_diff,NULL);
991:     PetscOptionsString("-rd_view_binary","File name to hold final solution","",rd->view_binary,rd->view_binary,sizeof(rd->view_binary),NULL);
992:   }
993:   PetscOptionsEnd();

995:   switch (rd->initial) {
996:   case 1:
997:   case 2:
998:     rd->rho            = 1.;
999:     rd->c              = 1.;
1000:     rd->K_R            = 1.;
1001:     rd->K_p            = 1.;
1002:     rd->sigma_b        = 0.25;
1003:     rd->MaterialEnergy = RDMaterialEnergy_Reduced;
1004:     break;
1005:   case 3:
1006:     /* Table 2 */
1007:     rd->rho     = 1.17e-3 * kilogram / (meter*meter*meter);                      /* density */
1008:     rd->K_R     = 7.44e18 * PetscPowRealInt(meter,5) * PetscPowReal(Kelvin,3.5) * PetscPowRealInt(kilogram,-2); /*  */
1009:     rd->K_p     = 2.33e20 * PetscPowRealInt(meter,5) * PetscPowReal(Kelvin,3.5) * PetscPowRealInt(kilogram,-2); /*  */
1010:     rd->I_H     = 2.179e-18 * Joule;                                             /* Hydrogen ionization potential */
1011:     rd->m_p     = 1.673e-27 * kilogram;                                          /* proton mass */
1012:     rd->m_e     = 9.109e-31 * kilogram;                                          /* electron mass */
1013:     rd->h       = 6.626e-34 * Joule * second;                                    /* Planck's constant */
1014:     rd->k       = 1.381e-23 * Joule / Kelvin;                                    /* Boltzman constant */
1015:     rd->c       = 3.00e8 * meter / second;                                       /* speed of light */
1016:     rd->sigma_b = 5.67e-8 * Watt * PetscPowRealInt(meter,-2) * PetscPowRealInt(Kelvin,-4);             /* Stefan-Boltzman constant */
1017:     rd->MaterialEnergy = RDMaterialEnergy_Saha;
1018:     break;
1019:   }

1021:   DMDACreate1d(comm,DM_BOUNDARY_NONE,20,sizeof(RDNode)/sizeof(PetscScalar),1,NULL,&rd->da);
1022:   DMSetFromOptions(rd->da);
1023:   DMSetUp(rd->da);
1024:   DMDASetFieldName(rd->da,0,"E");
1025:   DMDASetFieldName(rd->da,1,"T");
1026:   DMDASetUniformCoordinates(rd->da,0.,1.,0.,0.,0.,0.);

1028:   *inrd = rd;
1029:   return(0);
1030: }

1032: int main(int argc, char *argv[])
1033: {
1035:   RD             rd;
1036:   TS             ts;
1037:   SNES           snes;
1038:   Vec            X;
1039:   Mat            A,B;
1040:   PetscInt       steps;
1041:   PetscReal      ftime;

1043:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
1044:   RDCreate(PETSC_COMM_WORLD,&rd);
1045:   DMCreateGlobalVector(rd->da,&X);
1046:   DMSetMatType(rd->da,MATAIJ);
1047:   DMCreateMatrix(rd->da,&B);
1048:   RDInitialState(rd,X);

1050:   TSCreate(PETSC_COMM_WORLD,&ts);
1051:   TSSetProblemType(ts,TS_NONLINEAR);
1052:   TSSetType(ts,TSTHETA);
1053:   TSSetDM(ts,rd->da);
1054:   switch (rd->discretization) {
1055:   case DISCRETIZATION_FD:
1056:     TSSetIFunction(ts,NULL,RDIFunction_FD,rd);
1057:     if (rd->jacobian == JACOBIAN_ANALYTIC) TSSetIJacobian(ts,B,B,RDIJacobian_FD,rd);
1058:     break;
1059:   case DISCRETIZATION_FE:
1060:     TSSetIFunction(ts,NULL,RDIFunction_FE,rd);
1061:     if (rd->jacobian == JACOBIAN_ANALYTIC) TSSetIJacobian(ts,B,B,RDIJacobian_FE,rd);
1062:     break;
1063:   }
1064:   TSSetMaxTime(ts,rd->final_time);
1065:   TSSetTimeStep(ts,1e-3);
1066:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
1067:   TSSetFromOptions(ts);

1069:   A = B;
1070:   TSGetSNES(ts,&snes);
1071:   switch (rd->jacobian) {
1072:   case JACOBIAN_ANALYTIC:
1073:     break;
1074:   case JACOBIAN_MATRIXFREE:
1075:     break;
1076:   case JACOBIAN_FD_COLORING: {
1077:     SNESSetJacobian(snes,A,B,SNESComputeJacobianDefaultColor,0);
1078:   } break;
1079:   case JACOBIAN_FD_FULL:
1080:     SNESSetJacobian(snes,A,B,SNESComputeJacobianDefault,ts);
1081:     break;
1082:   }

1084:   if (rd->test_diff) {
1085:     RDTestDifferentiation(rd);
1086:   }
1087:   TSSolve(ts,X);
1088:   TSGetSolveTime(ts,&ftime);
1089:   TSGetStepNumber(ts,&steps);
1090:   PetscPrintf(PETSC_COMM_WORLD,"Steps %D  final time %g\n",steps,(double)ftime);
1091:   if (rd->view_draw) {
1092:     RDView(rd,X,PETSC_VIEWER_DRAW_WORLD);
1093:   }
1094:   if (rd->view_binary[0]) {
1095:     PetscViewer viewer;
1096:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,rd->view_binary,FILE_MODE_WRITE,&viewer);
1097:     RDView(rd,X,viewer);
1098:     PetscViewerDestroy(&viewer);
1099:   }
1100:   VecDestroy(&X);
1101:   MatDestroy(&B);
1102:   RDDestroy(&rd);
1103:   TSDestroy(&ts);
1104:   PetscFinalize();
1105:   return ierr;
1106: }
1107: /*TEST

1109:     test:
1110:       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
1111:       requires: !single

1113:     test:
1114:       suffix: 2
1115:       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
1116:       requires: !single

1118:     test:
1119:       suffix: 3
1120:       nsize: 2
1121:       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
1122:       requires: !single

1124: TEST*/