Actual source code: ex10.c
petsc-3.4.2 2013-07-02

      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) { rd->MaterialEnergy(rd,n,Em,dEm); }

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

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

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

    152: /* Eq 4 */
    153: static void RDDiffusionCoefficient(RD rd,PetscBool limit,RDNode *n,RDNode *nx,PetscScalar *D_R,RDNode *dD_R,RDNode *dxD_R)
    154: {
    155:   PetscScalar sigma_R,denom;
    156:   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;

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

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

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

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

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

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

    260: 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)
    261: {
    263:   PetscBool      istheta;

    266:   DMGetLocalVector(rd->da,X0loc);
    267:   DMGetLocalVector(rd->da,Xloc);
    268:   DMGetLocalVector(rd->da,Xloc_t);

    270:   DMGlobalToLocalBegin(rd->da,X,INSERT_VALUES,*Xloc);
    271:   DMGlobalToLocalEnd(rd->da,X,INSERT_VALUES,*Xloc);
    272:   DMGlobalToLocalBegin(rd->da,Xdot,INSERT_VALUES,*Xloc_t);
    273:   DMGlobalToLocalEnd(rd->da,Xdot,INSERT_VALUES,*Xloc_t);

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

    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)
    316: {
    318:   PetscInt       minloc;
    319:   PetscReal      min;

    322:   VecMin(X,&minloc,&min);
    323:   if (min < 0) {
    324:     SNES snes;
    325:     *in  = PETSC_FALSE;
    326:     TSGetSNES(ts,&snes);
    327:     SNESSetFunctionDomainError(snes);
    328:     PetscInfo3(ts,"Domain violation at %D field %D value %G\n",minloc/2,minloc%2,min);
    329:   } else *in = PETSC_TRUE;
    330:   return(0);
    331: }

    333: /* Energy and temperature must remain positive */
    334: #define RDCheckDomain(rd,ts,X) do {                                \
    335:     PetscErrorCode _ierr;                                          \
    336:     PetscBool      _in;                                                \
    337:     _RDCheckDomain_Private(rd,ts,X,&_in);CHKERRQ(_ierr);    \
    338:     if (!_in) return(0);                              \
    339:   } while (0)

    343: static PetscErrorCode RDIFunction_FD(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
    344: {
    346:   RD             rd = (RD)ctx;
    347:   RDNode         *x,*x0,*xdot,*f;
    348:   Vec            X0loc,Xloc,Xloc_t;
    349:   PetscReal      hx,Theta,dt;
    350:   DMDALocalInfo  info;
    351:   PetscInt       i;

    354:   RDCheckDomain(rd,ts,X);
    355:   RDGetLocalArrays(rd,ts,X,Xdot,&Theta,&dt,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
    356:   DMDAVecGetArray(rd->da,F,&f);
    357:   DMDAGetLocalInfo(rd->da,&info);
    358:   VecZeroEntries(F);

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

    362:   for (i=info.xs; i<info.xs+info.xm; i++) {
    363:     PetscReal   rho = rd->rho;
    364:     PetscScalar Em_t,rad;

    366:     rad = (1.-Theta)*RDRadiation(rd,&x0[i],0) + Theta*RDRadiation(rd,&x[i],0);
    367:     if (rd->endpoint) {
    368:       PetscScalar Em0,Em1;
    369:       RDMaterialEnergy(rd,&x0[i],&Em0,NULL);
    370:       RDMaterialEnergy(rd,&x[i],&Em1,NULL);
    371:       Em_t = (Em1 - Em0) / dt;
    372:     } else {
    373:       RDNode dEm;
    374:       RDMaterialEnergy(rd,&x[i],NULL,&dEm);
    375:       Em_t = dEm.E * xdot[i].E + dEm.T * xdot[i].T;
    376:     }
    377:     /* Residuals are multiplied by the volume element (hx).  */
    378:     /* The temperature equation does not have boundary conditions */
    379:     f[i].T = hx*(rho*Em_t + rad);

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

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

    414: static PetscErrorCode RDIJacobian_FD(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
    415: {
    417:   RD             rd = (RD)ctx;
    418:   RDNode         *x,*x0,*xdot;
    419:   Vec            X0loc,Xloc,Xloc_t;
    420:   PetscReal      hx,Theta,dt;
    421:   DMDALocalInfo  info;
    422:   PetscInt       i;

    425:   RDCheckDomain(rd,ts,X);
    426:   RDGetLocalArrays(rd,ts,X,Xdot,&Theta,&dt,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
    427:   DMDAGetLocalInfo(rd->da,&info);
    428:   hx   = rd->L / (info.mx-1);
    429:   MatZeroEntries(*B);

    431:   for (i=info.xs; i<info.xs+info.xm; i++) {
    432:     PetscInt                  col[3];
    433:     PetscReal                 rho = rd->rho;
    434:     PetscScalar /*Em_t,rad,*/ K[2][6];
    435:     RDNode                    dEm_t,drad;

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

    439:     if (rd->endpoint) {
    440:       PetscScalar Em0,Em1;
    441:       RDNode      dEm1;
    442:       RDMaterialEnergy(rd,&x0[i],&Em0,NULL);
    443:       RDMaterialEnergy(rd,&x[i],&Em1,&dEm1);
    444:       /*Em_t = (Em1 - Em0) / (Theta*dt);*/
    445:       dEm_t.E = dEm1.E / (Theta*dt);
    446:       dEm_t.T = dEm1.T / (Theta*dt);
    447:     } else {
    448:       const PetscScalar epsilon = x[i].T * PETSC_SQRT_MACHINE_EPSILON;
    449:       RDNode            n1;
    450:       RDNode            dEm,dEm1;
    451:       PetscScalar       Em_TT;

    453:       n1.E = x[i].E;
    454:       n1.T = x[i].T+epsilon;
    455:       RDMaterialEnergy(rd,&x[i],NULL,&dEm);
    456:       RDMaterialEnergy(rd,&n1,NULL,&dEm1);
    457:       /* The Jacobian needs another derivative.  We finite difference here instead of
    458:        * propagating second derivatives through the ionization model. */
    459:       Em_TT = (dEm1.T - dEm.T) / epsilon;
    460:       /*Em_t = dEm.E * xdot[i].E + dEm.T * xdot[i].T;*/
    461:       dEm_t.E = dEm.E * a;
    462:       dEm_t.T = dEm.T * a + Em_TT * xdot[i].T;
    463:     }

    465:     PetscMemzero(K,sizeof(K));
    466:     /* Residuals are multiplied by the volume element (hx).  */
    467:     if (i == 0) {
    468:       PetscScalar D,bcTheta = rd->bcmidpoint ? Theta : 1.;
    469:       RDNode      n, nx;
    470:       RDNode      dD,dxD;

    472:       n.E  = (1.-bcTheta)*x0[0].E + bcTheta*x[0].E;
    473:       n.T  = (1.-bcTheta)*x0[0].T + bcTheta*x[0].T;
    474:       nx.E = ((1.-bcTheta)*(x0[1].E-x0[0].E) + bcTheta*(x[1].E-x[0].E))/hx;
    475:       nx.T = ((1.-bcTheta)*(x0[1].T-x0[0].T) + bcTheta*(x[1].T-x[0].T))/hx;
    476:       switch (rd->leftbc) {
    477:       case BC_ROBIN:
    478:         RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D,&dD,&dxD);
    479:         K[0][1*2+0] = (bcTheta/Theta)*hx*(1. -2.*D*(-1./hx) - 2.*nx.E*dD.E + 2.*nx.E*dxD.E/hx);
    480:         K[0][1*2+1] = (bcTheta/Theta)*hx*(-2.*nx.E*dD.T);
    481:         K[0][2*2+0] = (bcTheta/Theta)*hx*(-2.*D*(1./hx) - 2.*nx.E*dD.E - 2.*nx.E*dxD.E/hx);
    482:         break;
    483:       case BC_NEUMANN:
    484:         K[0][1*2+0] = -1./Theta;
    485:         K[0][2*2+0] = 1./Theta;
    486:         break;
    487:       default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D",rd->initial);
    488:       }
    489:     } else if (i == info.mx-1) {
    490:       K[0][0*2+0] = -1./Theta;
    491:       K[0][1*2+0] = 1./Theta;
    492:     } else {
    493:       /*PetscScalar diff;*/
    494:       RDNode ddiff[3];
    495:       /*diff = (1.-Theta)*RDDiffusion(rd,hx,x0,i,0) + Theta* */ RDDiffusion(rd,hx,x,i,ddiff);
    496:       K[0][0*2+0] = -hx*ddiff[0].E;
    497:       K[0][0*2+1] = -hx*ddiff[0].T;
    498:       K[0][1*2+0] = hx*(a - ddiff[1].E - drad.E);
    499:       K[0][1*2+1] = hx*(-ddiff[1].T - drad.T);
    500:       K[0][2*2+0] = -hx*ddiff[2].E;
    501:       K[0][2*2+1] = -hx*ddiff[2].T;
    502:     }

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

    507:     col[0] = i-1;
    508:     col[1] = i;
    509:     col[2] = i+1<info.mx ? i+1 : -1;
    510:     MatSetValuesBlocked(*B,1,&i,3,col,&K[0][0],INSERT_VALUES);
    511:   }
    512:   RDRestoreLocalArrays(rd,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
    513:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
    514:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
    515:   if (*A != *B) {
    516:     MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
    517:     MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
    518:   }
    519:   *mstr = SAME_NONZERO_PATTERN;
    520:   return(0);
    521: }

    523: /* Evaluate interpolants and derivatives at a select quadrature point */
    524: static void RDEvaluate(PetscReal interp[][2],PetscReal deriv[][2],PetscInt q,const RDNode x[],PetscInt i,RDNode *n,RDNode *nx)
    525: {
    526:   PetscInt j;
    527:   n->E = 0; n->T = 0; nx->E = 0; nx->T = 0;
    528:   for (j=0; j<2; j++) {
    529:     n->E  += interp[q][j] * x[i+j].E;
    530:     n->T  += interp[q][j] * x[i+j].T;
    531:     nx->E += deriv[q][j] * x[i+j].E;
    532:     nx->T += deriv[q][j] * x[i+j].T;
    533:   }
    534: }

    538: /*
    539:  Various quadrature rules.  The nonlinear terms are non-polynomial so no standard quadrature will be exact.
    540: */
    541: static PetscErrorCode RDGetQuadrature(RD rd,PetscReal hx,PetscInt *nq,PetscReal weight[],PetscReal interp[][2],PetscReal deriv[][2])
    542: {
    543:   PetscInt        q,j;
    544:   const PetscReal *refweight,(*refinterp)[2],(*refderiv)[2];

    547:   switch (rd->quadrature) {
    548:   case QUADRATURE_GAUSS1: {
    549:     static const PetscReal ww[1] = {1.},ii[1][2] = {{0.5,0.5}},dd[1][2] = {{-1.,1.}};
    550:     *nq = 1; refweight = ww; refinterp = ii; refderiv = dd;
    551:   } break;
    552:   case QUADRATURE_GAUSS2: {
    553:     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};
    554:     *nq = 2; refweight = ww; refinterp = ii; refderiv = dd;
    555:   } break;
    556:   case QUADRATURE_GAUSS3: {
    557:     static const PetscReal ii[3][2] = {{0.8872983346207417,0.1127016653792583},{0.5,0.5},{0.1127016653792583,0.8872983346207417}},
    558:                            dd[3][2] = {{-1,1},{-1,1},{-1,1}},ww[3] = {5./18,8./18,5./18};
    559:     *nq = 3; refweight = ww; refinterp = ii; refderiv = dd;
    560:   } break;
    561:   case QUADRATURE_GAUSS4: {
    562:     static const PetscReal ii[][2] = {{0.93056815579702623,0.069431844202973658},
    563:                                       {0.66999052179242813,0.33000947820757187},
    564:                                       {0.33000947820757187,0.66999052179242813},
    565:                                       {0.069431844202973658,0.93056815579702623}},
    566:                            dd[][2] = {{-1,1},{-1,1},{-1,1},{-1,1}},ww[] = {0.17392742256872692,0.3260725774312731,0.3260725774312731,0.17392742256872692};

    568:     *nq = 4; refweight = ww; refinterp = ii; refderiv = dd;
    569:   } break;
    570:   case QUADRATURE_LOBATTO2: {
    571:     static const PetscReal ii[2][2] = {{1.,0.},{0.,1.}},dd[2][2] = {{-1.,1.},{-1.,1.}},ww[2] = {0.5,0.5};
    572:     *nq = 2; refweight = ww; refinterp = ii; refderiv = dd;
    573:   } break;
    574:   case QUADRATURE_LOBATTO3: {
    575:     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};
    576:     *nq = 3; refweight = ww; refinterp = ii; refderiv = dd;
    577:   } break;
    578:   default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown quadrature %d",(int)rd->quadrature);
    579:   }

    581:   for (q=0; q<*nq; q++) {
    582:     weight[q] = refweight[q] * hx;
    583:     for (j=0; j<2; j++) {
    584:       interp[q][j] = refinterp[q][j];
    585:       deriv[q][j]  = refderiv[q][j] / hx;
    586:     }
    587:   }
    588:   return(0);
    589: }

    593: /*
    594:  Finite element version
    595: */
    596: static PetscErrorCode RDIFunction_FE(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
    597: {
    599:   RD             rd = (RD)ctx;
    600:   RDNode         *x,*x0,*xdot,*f;
    601:   Vec            X0loc,Xloc,Xloc_t,Floc;
    602:   PetscReal      hx,Theta,dt,weight[5],interp[5][2],deriv[5][2];
    603:   DMDALocalInfo  info;
    604:   PetscInt       i,j,q,nq;

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

    610:   DMGetLocalVector(rd->da,&Floc);
    611:   VecZeroEntries(Floc);
    612:   DMDAVecGetArray(rd->da,Floc,&f);
    613:   DMDAGetLocalInfo(rd->da,&info);

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

    619:   for (i=info.xs; i<PetscMin(info.xs+info.xm,info.mx-1); i++) {
    620:     for (q=0; q<nq; q++) {
    621:       PetscReal   rho = rd->rho;
    622:       PetscScalar Em_t,rad,D_R,D0_R;
    623:       RDNode      n,n0,nx,n0x,nt,ntx;
    624:       RDEvaluate(interp,deriv,q,x,i,&n,&nx);
    625:       RDEvaluate(interp,deriv,q,x0,i,&n0,&n0x);
    626:       RDEvaluate(interp,deriv,q,xdot,i,&nt,&ntx);

    628:       rad = (1.-Theta)*RDRadiation(rd,&n0,0) + Theta*RDRadiation(rd,&n,0);
    629:       if (rd->endpoint) {
    630:         PetscScalar Em0,Em1;
    631:         RDMaterialEnergy(rd,&n0,&Em0,NULL);
    632:         RDMaterialEnergy(rd,&n,&Em1,NULL);
    633:         Em_t = (Em1 - Em0) / dt;
    634:       } else {
    635:         RDNode dEm;
    636:         RDMaterialEnergy(rd,&n,NULL,&dEm);
    637:         Em_t = dEm.E * nt.E + dEm.T * nt.T;
    638:       }
    639:       RDDiffusionCoefficient(rd,PETSC_TRUE,&n0,&n0x,&D0_R,0,0);
    640:       RDDiffusionCoefficient(rd,PETSC_TRUE,&n,&nx,&D_R,0,0);
    641:       for (j=0; j<2; j++) {
    642:         f[i+j].E += (deriv[q][j] * weight[q] * ((1.-Theta)*D0_R*n0x.E + Theta*D_R*nx.E)
    643:                      + interp[q][j] * weight[q] * (nt.E - rad));
    644:         f[i+j].T += interp[q][j] * weight[q] * (rho * Em_t + rad);
    645:       }
    646:     }
    647:   }
    648:   if (info.xs == 0) {
    649:     switch (rd->leftbc) {
    650:     case BC_ROBIN: {
    651:       PetscScalar D_R,D_R_bc;
    652:       PetscReal   ratio,bcTheta = rd->bcmidpoint ? Theta : 1.;
    653:       RDNode      n, nx;

    655:       n.E  = (1-bcTheta)*x0[0].E + bcTheta*x[0].E;
    656:       n.T  = (1-bcTheta)*x0[0].T + bcTheta*x[0].T;
    657:       nx.E = (x[1].E-x[0].E)/hx;
    658:       nx.T = (x[1].T-x[0].T)/hx;
    659:       RDDiffusionCoefficient(rd,PETSC_TRUE,&n,&nx,&D_R,0,0);
    660:       RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D_R_bc,0,0);
    661:       ratio = PetscRealPart(D_R/D_R_bc);
    662:       if (ratio > 1.) SETERRQ(PETSC_COMM_SELF,1,"Limited diffusivity is greater than unlimited");
    663:       if (ratio < 1e-3) SETERRQ(PETSC_COMM_SELF,1,"Heavily limited diffusivity");
    664:       f[0].E += -ratio*0.5*(rd->Eapplied - n.E);
    665:     } break;
    666:     case BC_NEUMANN:
    667:       /* homogeneous Neumann is the natural condition */
    668:       break;
    669:     default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D",rd->initial);
    670:     }
    671:   }

    673:   RDRestoreLocalArrays(rd,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
    674:   DMDAVecRestoreArray(rd->da,Floc,&f);
    675:   VecZeroEntries(F);
    676:   DMLocalToGlobalBegin(rd->da,Floc,ADD_VALUES,F);
    677:   DMLocalToGlobalEnd(rd->da,Floc,ADD_VALUES,F);
    678:   DMRestoreLocalVector(rd->da,&Floc);

    680:   if (rd->monitor_residual) {RDStateView(rd,X,Xdot,F);}
    681:   return(0);
    682: }

    686: static PetscErrorCode RDIJacobian_FE(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
    687: {
    689:   RD             rd = (RD)ctx;
    690:   RDNode         *x,*x0,*xdot;
    691:   Vec            X0loc,Xloc,Xloc_t;
    692:   PetscReal      hx,Theta,dt,weight[5],interp[5][2],deriv[5][2];
    693:   DMDALocalInfo  info;
    694:   PetscInt       i,j,k,q,nq;
    695:   PetscScalar    K[4][4];

    698:   RDCheckDomain(rd,ts,X);
    699:   RDGetLocalArrays(rd,ts,X,Xdot,&Theta,&dt,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
    700:   DMDAGetLocalInfo(rd->da,&info);
    701:   hx   = rd->L / (info.mx-1);
    702:   RDGetQuadrature(rd,hx,&nq,weight,interp,deriv);
    703:   MatZeroEntries(*B);
    704:   for (i=info.xs; i<PetscMin(info.xs+info.xm,info.mx-1); i++) {
    705:     PetscInt rc[2];

    707:     rc[0] = i; rc[1] = i+1;
    708:     PetscMemzero(K,sizeof(K));
    709:     for (q=0; q<nq; q++) {
    710:       PetscScalar D_R,PETSC_UNUSED rad;
    711:       RDNode      n,nx,nt,ntx,drad,dD_R,dxD_R,dEm;
    712:       RDEvaluate(interp,deriv,q,x,i,&n,&nx);
    713:       RDEvaluate(interp,deriv,q,xdot,i,&nt,&ntx);
    714:       rad = RDRadiation(rd,&n,&drad);
    715:       RDDiffusionCoefficient(rd,PETSC_TRUE,&n,&nx,&D_R,&dD_R,&dxD_R);
    716:       RDMaterialEnergy(rd,&n,NULL,&dEm);
    717:       for (j=0; j<2; j++) {
    718:         for (k=0; k<2; k++) {
    719:           K[j*2+0][k*2+0] += (+interp[q][j] * weight[q] * (a - drad.E) * interp[q][k]
    720:                               + deriv[q][j] * weight[q] * ((D_R + dxD_R.E * nx.E) * deriv[q][k] + dD_R.E * nx.E * interp[q][k]));
    721:           K[j*2+0][k*2+1] += (+interp[q][j] * weight[q] * (-drad.T * interp[q][k])
    722:                               + deriv[q][j] * weight[q] * (dxD_R.T * deriv[q][k] + dD_R.T * interp[q][k]) * nx.E);
    723:           K[j*2+1][k*2+0] +=   interp[q][j] * weight[q] * drad.E * interp[q][k];
    724:           K[j*2+1][k*2+1] +=   interp[q][j] * weight[q] * (a * rd->rho * dEm.T + drad.T) * interp[q][k];
    725:         }
    726:       }
    727:     }
    728:     MatSetValuesBlocked(*B,2,rc,2,rc,&K[0][0],ADD_VALUES);
    729:   }
    730:   if (info.xs == 0) {
    731:     switch (rd->leftbc) {
    732:     case BC_ROBIN: {
    733:       PetscScalar D_R,D_R_bc;
    734:       PetscReal   ratio;
    735:       RDNode      n, nx;

    737:       n.E  = (1-Theta)*x0[0].E + Theta*x[0].E;
    738:       n.T  = (1-Theta)*x0[0].T + Theta*x[0].T;
    739:       nx.E = (x[1].E-x[0].E)/hx;
    740:       nx.T = (x[1].T-x[0].T)/hx;
    741:       RDDiffusionCoefficient(rd,PETSC_TRUE,&n,&nx,&D_R,0,0);
    742:       RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D_R_bc,0,0);
    743:       ratio = PetscRealPart(D_R/D_R_bc);
    744:       MatSetValue(*B,0,0,ratio*0.5,ADD_VALUES);
    745:     } break;
    746:     case BC_NEUMANN:
    747:       /* homogeneous Neumann is the natural condition */
    748:       break;
    749:     default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %D",rd->initial);
    750:     }
    751:   }

    753:   RDRestoreLocalArrays(rd,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);
    754:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
    755:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
    756:   if (*A != *B) {
    757:     MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
    758:     MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
    759:   }
    760:   *mstr = SAME_NONZERO_PATTERN;
    761:   return(0);
    762: }

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

    769: static PetscErrorCode RDInitialState(RD rd,Vec X)
    770: {
    771:   DMDALocalInfo  info;
    772:   PetscInt       i;
    773:   RDNode         *x;

    777:   DMDAGetLocalInfo(rd->da,&info);
    778:   DMDAVecGetArray(rd->da,X,&x);
    779:   for (i=info.xs; i<info.xs+info.xm; i++) {
    780:     PetscReal coord = i*rd->L/(info.mx-1);
    781:     switch (rd->initial) {
    782:     case 1:
    783:       x[i].E = 0.001;
    784:       x[i].T = RDRadiationTemperature(rd,x[i].E);
    785:       break;
    786:     case 2:
    787:       x[i].E = 0.001 + 100.*PetscExpScalar(-PetscSqr(coord/0.1));
    788:       x[i].T = RDRadiationTemperature(rd,x[i].E);
    789:       break;
    790:     case 3:
    791:       x[i].E = 7.56e-2 * rd->unit.Joule / pow(rd->unit.meter,3);
    792:       x[i].T = RDRadiationTemperature(rd,x[i].E);
    793:       break;
    794:     default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No initial state %D",rd->initial);
    795:     }
    796:   }
    797:   DMDAVecRestoreArray(rd->da,X,&x);
    798:   return(0);
    799: }

    803: static PetscErrorCode RDView(RD rd,Vec X,PetscViewer viewer)
    804: {
    806:   Vec            Y;
    807:   RDNode         *x;
    808:   PetscScalar    *y;
    809:   PetscInt       i,m,M;
    810:   const PetscInt *lx;
    811:   DM             da;
    812:   MPI_Comm       comm;

    815:   /*
    816:     Create a DMDA (one dof per node, zero stencil width, same layout) to hold Trad
    817:     (radiation temperature).  It is not necessary to create a DMDA for this, but this way
    818:     output and visualization will have meaningful variable names and correct scales.
    819:   */
    820:   DMDAGetInfo(rd->da,0, &M,0,0, 0,0,0, 0,0,0,0,0,0);
    821:   DMDAGetOwnershipRanges(rd->da,&lx,0,0);
    822:   PetscObjectGetComm((PetscObject)rd->da,&comm);
    823:   DMDACreate1d(comm,DMDA_BOUNDARY_NONE,M,1,0,lx,&da);
    824:   DMDASetUniformCoordinates(da,0.,rd->L,0.,0.,0.,0.);
    825:   DMDASetFieldName(da,0,"T_rad");
    826:   DMCreateGlobalVector(da,&Y);

    828:   /* Compute the radiation temperature from the solution at each node */
    829:   VecGetLocalSize(Y,&m);
    830:   VecGetArray(X,(PetscScalar**)&x);
    831:   VecGetArray(Y,&y);
    832:   for (i=0; i<m; i++) y[i] = RDRadiationTemperature(rd,x[i].E);
    833:   VecRestoreArray(X,(PetscScalar**)&x);
    834:   VecRestoreArray(Y,&y);

    836:   VecView(Y,viewer);
    837:   VecDestroy(&Y);
    838:   DMDestroy(&da);
    839:   return(0);
    840: }

    844: static PetscErrorCode RDTestDifferentiation(RD rd)
    845: {
    846:   MPI_Comm       comm;
    848:   RDNode         n,nx;
    849:   PetscScalar    epsilon;

    852:   PetscObjectGetComm((PetscObject)rd->da,&comm);
    853:   epsilon = 1e-8;
    854:   {
    855:     RDNode      dEm,fdEm;
    856:     PetscScalar T0 = 1000.,T1 = T0*(1.+epsilon),Em0,Em1;
    857:     n.E = 1.;
    858:     n.T = T0;
    859:     rd->MaterialEnergy(rd,&n,&Em0,&dEm);
    860:     n.E = 1.+epsilon;
    861:     n.T = T0;
    862:     rd->MaterialEnergy(rd,&n,&Em1,0);
    863:     fdEm.E = (Em1-Em0)/epsilon;
    864:     n.E = 1.;
    865:     n.T = T1;
    866:     rd->MaterialEnergy(rd,&n,&Em1,0);
    867:     fdEm.T = (Em1-Em0)/(T0*epsilon);
    868:     PetscPrintf(comm,"dEm {%G,%G}, fdEm {%G,%G}, diff {%G,%G}\n",PetscRealPart(dEm.E),PetscRealPart(dEm.T),
    869:                        PetscRealPart(fdEm.E),PetscRealPart(fdEm.T),PetscRealPart(dEm.E-fdEm.E),PetscRealPart(dEm.T-fdEm.T));
    870:   }
    871:   {
    872:     PetscScalar D0,D;
    873:     RDNode      dD,dxD,fdD,fdxD;
    874:     n.E = 1.;          n.T = 1.;           nx.E = 1.;          n.T = 1.;
    875:     RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D0,&dD,&dxD);
    876:     n.E = 1.+epsilon;  n.T = 1.;           nx.E = 1.;          n.T = 1.;
    877:     RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D,0,0); fdD.E = (D-D0)/epsilon;
    878:     n.E = 1;           n.T = 1.+epsilon;   nx.E = 1.;          n.T = 1.;
    879:     RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D,0,0); fdD.T = (D-D0)/epsilon;
    880:     n.E = 1;           n.T = 1.;           nx.E = 1.+epsilon;  n.T = 1.;
    881:     RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D,0,0); fdxD.E = (D-D0)/epsilon;
    882:     n.E = 1;           n.T = 1.;           nx.E = 1.;          n.T = 1.+epsilon;
    883:     RDDiffusionCoefficient(rd,rd->bclimit,&n,&nx,&D,0,0); fdxD.T = (D-D0)/epsilon;
    884:     PetscPrintf(comm,"dD {%G,%G}, fdD {%G,%G}, diff {%G,%G}\n",PetscRealPart(dD.E),PetscRealPart(dD.T),
    885:                        PetscRealPart(fdD.E),PetscRealPart(fdD.T),PetscRealPart(dD.E-fdD.E),PetscRealPart(dD.T-fdD.T));
    886:     PetscPrintf(comm,"dxD {%G,%G}, fdxD {%G,%G}, diffx {%G,%G}\n",PetscRealPart(dxD.E),PetscRealPart(dxD.T),
    887:                        PetscRealPart(fdxD.E),PetscRealPart(fdxD.T),PetscRealPart(dxD.E-fdxD.E),PetscRealPart(dxD.T-fdxD.T));
    888:   }
    889:   {
    890:     PetscInt    i;
    891:     PetscReal   hx = 1.;
    892:     PetscScalar a0;
    893:     RDNode      n0[3] = {{1.,1.},{5.,3.},{4.,2.}},n1[3],d[3],fd[3];
    894:     a0 = RDDiffusion(rd,hx,n0,1,d);
    895:     for (i=0; i<3; i++) {
    896:       PetscMemcpy(n1,n0,sizeof(n0)); n1[i].E += epsilon;
    897:       fd[i].E = (RDDiffusion(rd,hx,n1,1,0)-a0)/epsilon;
    898:       PetscMemcpy(n1,n0,sizeof(n0)); n1[i].T += epsilon;
    899:       fd[i].T = (RDDiffusion(rd,hx,n1,1,0)-a0)/epsilon;
    900:       PetscPrintf(comm,"ddiff[%D] {%G,%G}, fd {%G %G}, diff {%G,%G}\n",i,PetscRealPart(d[i].E),PetscRealPart(d[i].T),
    901:                             PetscRealPart(fd[i].E),PetscRealPart(fd[i].T),PetscRealPart(d[i].E-fd[i].E),PetscRealPart(d[i].T-fd[i].T));
    902:     }
    903:   }
    904:   {
    905:     PetscScalar rad0,rad;
    906:     RDNode      drad,fdrad;
    907:     n.E  = 1.;         n.T = 1.;
    908:     rad0 = RDRadiation(rd,&n,&drad);
    909:     n.E  = 1.+epsilon; n.T = 1.;
    910:     rad  = RDRadiation(rd,&n,0); fdrad.E = (rad-rad0)/epsilon;
    911:     n.E  = 1.;         n.T = 1.+epsilon;
    912:     rad  = RDRadiation(rd,&n,0); fdrad.T = (rad-rad0)/epsilon;
    913:     PetscPrintf(comm,"drad {%G,%G}, fdrad {%G,%G}, diff {%G,%G}\n",PetscRealPart(drad.E),PetscRealPart(drad.T),
    914:                        PetscRealPart(fdrad.E),PetscRealPart(fdrad.T),PetscRealPart(drad.E-drad.E),PetscRealPart(drad.T-fdrad.T));
    915:   }
    916:   return(0);
    917: }

    921: static PetscErrorCode RDCreate(MPI_Comm comm,RD *inrd)
    922: {
    924:   RD             rd;
    925:   PetscReal      meter,kilogram,second,Kelvin,Joule,Watt;

    928:   *inrd = 0;
    929:   PetscNew(struct _n_RD,&rd);

    931:   PetscOptionsBegin(comm,NULL,"Options for nonequilibrium radiation-diffusion with RD ionization",NULL);
    932:   {
    933:     rd->initial = 1;
    934:     PetscOptionsInt("-rd_initial","Initial condition (1=Marshak, 2=Blast, 3=Marshak+)","",rd->initial,&rd->initial,0);
    935:     switch (rd->initial) {
    936:     case 1:
    937:     case 2:
    938:       rd->unit.kilogram = 1.;
    939:       rd->unit.meter    = 1.;
    940:       rd->unit.second   = 1.;
    941:       rd->unit.Kelvin   = 1.;
    942:       break;
    943:     case 3:
    944:       rd->unit.kilogram = 1.e12;
    945:       rd->unit.meter    = 1.;
    946:       rd->unit.second   = 1.e9;
    947:       rd->unit.Kelvin   = 1.;
    948:       break;
    949:     default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown initial condition %d",rd->initial);
    950:     }
    951:     /* Fundamental units */
    952:     PetscOptionsReal("-rd_unit_meter","Length of 1 meter in nondimensional units","",rd->unit.meter,&rd->unit.meter,0);
    953:     PetscOptionsReal("-rd_unit_kilogram","Mass of 1 kilogram in nondimensional units","",rd->unit.kilogram,&rd->unit.kilogram,0);
    954:     PetscOptionsReal("-rd_unit_second","Time of a second in nondimensional units","",rd->unit.second,&rd->unit.second,0);
    955:     PetscOptionsReal("-rd_unit_Kelvin","Temperature of a Kelvin in nondimensional units","",rd->unit.Kelvin,&rd->unit.Kelvin,0);
    956:     /* Derived units */
    957:     rd->unit.Joule = rd->unit.kilogram*PetscSqr(rd->unit.meter/rd->unit.second);
    958:     rd->unit.Watt  = rd->unit.Joule/rd->unit.second;
    959:     /* Local aliases */
    960:     meter    = rd->unit.meter;
    961:     kilogram = rd->unit.kilogram;
    962:     second   = rd->unit.second;
    963:     Kelvin   = rd->unit.Kelvin;
    964:     Joule    = rd->unit.Joule;
    965:     Watt     = rd->unit.Watt;

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

    1014:   switch (rd->initial) {
    1015:   case 1:
    1016:   case 2:
    1017:     rd->rho            = 1.;
    1018:     rd->c              = 1.;
    1019:     rd->K_R            = 1.;
    1020:     rd->K_p            = 1.;
    1021:     rd->sigma_b        = 0.25;
    1022:     rd->MaterialEnergy = RDMaterialEnergy_Reduced;
    1023:     break;
    1024:   case 3:
    1025:     /* Table 2 */
    1026:     rd->rho     = 1.17e-3 * kilogram / (meter*meter*meter);                      /* density */
    1027:     rd->K_R     = 7.44e18 * pow(meter,5.) * pow(Kelvin,3.5) * pow(kilogram,-2.); /*  */
    1028:     rd->K_p     = 2.33e20 * pow(meter,5.) * pow(Kelvin,3.5) * pow(kilogram,-2.); /*  */
    1029:     rd->I_H     = 2.179e-18 * Joule;                                             /* Hydrogen ionization potential */
    1030:     rd->m_p     = 1.673e-27 * kilogram;                                          /* proton mass */
    1031:     rd->m_e     = 9.109e-31 * kilogram;                                          /* electron mass */
    1032:     rd->h       = 6.626e-34 * Joule * second;                                    /* Planck's constant */
    1033:     rd->k       = 1.381e-23 * Joule / Kelvin;                                    /* Boltzman constant */
    1034:     rd->c       = 3.00e8 * meter / second;                                       /* speed of light */
    1035:     rd->sigma_b = 5.67e-8 * Watt * pow(meter,-2.) * pow(Kelvin,-4.);             /* Stefan-Boltzman constant */
    1036:     rd->MaterialEnergy = RDMaterialEnergy_Saha;
    1037:     break;
    1038:   }

    1040:   DMDACreate1d(comm,DMDA_BOUNDARY_NONE,-20,sizeof(RDNode)/sizeof(PetscScalar),1,NULL,&rd->da);
    1041:   DMDASetFieldName(rd->da,0,"E");
    1042:   DMDASetFieldName(rd->da,1,"T");
    1043:   DMDASetUniformCoordinates(rd->da,0.,1.,0.,0.,0.,0.);

    1045:   *inrd = rd;
    1046:   return(0);
    1047: }

    1051: int main(int argc, char *argv[])
    1052: {
    1054:   RD             rd;
    1055:   TS             ts;
    1056:   SNES           snes;
    1057:   Vec            X;
    1058:   Mat            A,B;
    1059:   PetscInt       steps;
    1060:   PetscReal      ftime;

    1062:   PetscInitialize(&argc,&argv,0,help);
    1063:   RDCreate(PETSC_COMM_WORLD,&rd);
    1064:   DMCreateGlobalVector(rd->da,&X);
    1065:   DMCreateMatrix(rd->da,MATAIJ,&B);
    1066:   RDInitialState(rd,X);

    1068:   TSCreate(PETSC_COMM_WORLD,&ts);
    1069:   TSSetProblemType(ts,TS_NONLINEAR);
    1070:   TSSetType(ts,TSTHETA);
    1071:   TSSetDM(ts,rd->da);
    1072:   switch (rd->discretization) {
    1073:   case DISCRETIZATION_FD:
    1074:     TSSetIFunction(ts,NULL,RDIFunction_FD,rd);
    1075:     if (rd->jacobian == JACOBIAN_ANALYTIC) TSSetIJacobian(ts,B,B,RDIJacobian_FD,rd);
    1076:     break;
    1077:   case DISCRETIZATION_FE:
    1078:     TSSetIFunction(ts,NULL,RDIFunction_FE,rd);
    1079:     if (rd->jacobian == JACOBIAN_ANALYTIC) TSSetIJacobian(ts,B,B,RDIJacobian_FE,rd);
    1080:     break;
    1081:   }
    1082:   TSSetDuration(ts,10000,rd->final_time);
    1083:   TSSetInitialTimeStep(ts,0.,1e-3);
    1084:   TSSetFromOptions(ts);

    1086:   A = B;
    1087:   TSGetSNES(ts,&snes);
    1088:   switch (rd->jacobian) {
    1089:   case JACOBIAN_ANALYTIC:
    1090:     break;
    1091:   case JACOBIAN_MATRIXFREE:
    1092:     break;
    1093:   case JACOBIAN_FD_COLORING: {
    1094:     SNESSetJacobian(snes,A,B,SNESComputeJacobianDefaultColor,0);
    1095:   } break;
    1096:   case JACOBIAN_FD_FULL:
    1097:     SNESSetJacobian(snes,A,B,SNESComputeJacobianDefault,ts);
    1098:     break;
    1099:   }

    1101:   if (rd->test_diff) {
    1102:     RDTestDifferentiation(rd);
    1103:   }
    1104:   TSSolve(ts,X);
    1105:   TSGetSolveTime(ts,&ftime);
    1106:   TSGetTimeStepNumber(ts,&steps);
    1107:   PetscPrintf(PETSC_COMM_WORLD,"Steps %D  final time %G\n",steps,ftime);
    1108:   if (rd->view_draw) {
    1109:     RDView(rd,X,PETSC_VIEWER_DRAW_WORLD);
    1110:   }
    1111:   if (rd->view_binary[0]) {
    1112:     PetscViewer viewer;
    1113:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,rd->view_binary,FILE_MODE_WRITE,&viewer);
    1114:     RDView(rd,X,viewer);
    1115:     PetscViewerDestroy(&viewer);
    1116:   }
    1117:   VecDestroy(&X);
    1118:   MatDestroy(&B);
    1119:   RDDestroy(&rd);
    1120:   TSDestroy(&ts);
    1121:   PetscFinalize();
    1122:   return 0;
    1123: }