1 static const char help[] = "1D nonequilibrium radiation diffusion with Saha ionization model.\n\n"; 2 3 /* 4 This example implements the model described in 5 6 Rauenzahn, Mousseau, Knoll. "Temporal accuracy of the nonequilibrium radiation diffusion 7 equations employing a Saha ionization model" 2005. 8 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 */ 13 14 #include <petscts.h> 15 #include <petscdm.h> 16 #include <petscdmda.h> 17 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}; 26 27 typedef struct { 28 PetscScalar E; /* radiation energy */ 29 PetscScalar T; /* material temperature */ 30 } RDNode; 31 32 typedef struct { 33 PetscReal meter,kilogram,second,Kelvin; /* Fundamental units */ 34 PetscReal Joule,Watt; /* Derived units */ 35 } RDUnit; 36 37 typedef struct _n_RD *RD; 38 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; 55 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; 58 59 /* Domain and boundary conditions */ 60 PetscReal Eapplied; /* Radiation flux from the left */ 61 PetscReal L; /* Length of domain */ 62 PetscReal final_time; 63 }; 64 65 static PetscErrorCode RDDestroy(RD *rd) 66 { 67 PetscErrorCode ierr; 68 69 PetscFunctionBeginUser; 70 ierr = DMDestroy(&(*rd)->da);CHKERRQ(ierr); 71 ierr = PetscFree(*rd);CHKERRQ(ierr); 72 PetscFunctionReturn(0); 73 } 74 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); } 84 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 } 98 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 } 142 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 } 150 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; 156 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 } 173 174 static PetscErrorCode RDStateView(RD rd,Vec X,Vec Xdot,Vec F) 175 { 176 PetscErrorCode ierr; 177 DMDALocalInfo info; 178 PetscInt i; 179 const RDNode *x,*xdot,*f; 180 MPI_Comm comm; 181 182 PetscFunctionBeginUser; 183 ierr = PetscObjectGetComm((PetscObject)rd->da,&comm);CHKERRQ(ierr); 184 ierr = DMDAGetLocalInfo(rd->da,&info);CHKERRQ(ierr); 185 ierr = DMDAVecGetArrayRead(rd->da,X,(void*)&x);CHKERRQ(ierr); 186 ierr = DMDAVecGetArrayRead(rd->da,Xdot,(void*)&xdot);CHKERRQ(ierr); 187 ierr = DMDAVecGetArrayRead(rd->da,F,(void*)&f);CHKERRQ(ierr); 188 for (i=info.xs; i<info.xs+info.xm; i++) { 189 ierr = 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));CHKERRQ(ierr); 191 } 192 ierr = DMDAVecRestoreArrayRead(rd->da,X,(void*)&x);CHKERRQ(ierr); 193 ierr = DMDAVecRestoreArrayRead(rd->da,Xdot,(void*)&xdot);CHKERRQ(ierr); 194 ierr = DMDAVecRestoreArrayRead(rd->da,F,(void*)&f);CHKERRQ(ierr); 195 ierr = PetscSynchronizedFlush(comm,PETSC_STDOUT);CHKERRQ(ierr); 196 PetscFunctionReturn(0); 197 } 198 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 } 215 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; 221 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; 232 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; 243 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 } 254 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 { 257 PetscErrorCode ierr; 258 PetscBool istheta; 259 260 PetscFunctionBeginUser; 261 ierr = DMGetLocalVector(rd->da,X0loc);CHKERRQ(ierr); 262 ierr = DMGetLocalVector(rd->da,Xloc);CHKERRQ(ierr); 263 ierr = DMGetLocalVector(rd->da,Xloc_t);CHKERRQ(ierr); 264 265 ierr = DMGlobalToLocalBegin(rd->da,X,INSERT_VALUES,*Xloc);CHKERRQ(ierr); 266 ierr = DMGlobalToLocalEnd(rd->da,X,INSERT_VALUES,*Xloc);CHKERRQ(ierr); 267 ierr = DMGlobalToLocalBegin(rd->da,Xdot,INSERT_VALUES,*Xloc_t);CHKERRQ(ierr); 268 ierr = DMGlobalToLocalEnd(rd->da,Xdot,INSERT_VALUES,*Xloc_t);CHKERRQ(ierr); 269 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 ierr = PetscObjectTypeCompare((PetscObject)ts,TSTHETA,&istheta);CHKERRQ(ierr); 276 if (istheta && rd->endpoint) { 277 ierr = TSThetaGetTheta(ts,Theta);CHKERRQ(ierr); 278 } else *Theta = 1.; 279 280 ierr = TSGetTimeStep(ts,dt);CHKERRQ(ierr); 281 ierr = VecWAXPY(*X0loc,-(*Theta)*(*dt),*Xloc_t,*Xloc);CHKERRQ(ierr); /* back out the value at the start of this step */ 282 if (rd->endpoint) { 283 ierr = VecWAXPY(*Xloc,*dt,*Xloc_t,*X0loc);CHKERRQ(ierr); /* move the abscissa to the end of the step */ 284 } 285 286 ierr = DMDAVecGetArray(rd->da,*X0loc,x0);CHKERRQ(ierr); 287 ierr = DMDAVecGetArray(rd->da,*Xloc,x);CHKERRQ(ierr); 288 ierr = DMDAVecGetArray(rd->da,*Xloc_t,xdot);CHKERRQ(ierr); 289 PetscFunctionReturn(0); 290 } 291 292 static PetscErrorCode RDRestoreLocalArrays(RD rd,Vec *X0loc,RDNode **x0,Vec *Xloc,RDNode **x,Vec *Xloc_t,RDNode **xdot) 293 { 294 PetscErrorCode ierr; 295 296 PetscFunctionBeginUser; 297 ierr = DMDAVecRestoreArray(rd->da,*X0loc,x0);CHKERRQ(ierr); 298 ierr = DMDAVecRestoreArray(rd->da,*Xloc,x);CHKERRQ(ierr); 299 ierr = DMDAVecRestoreArray(rd->da,*Xloc_t,xdot);CHKERRQ(ierr); 300 ierr = DMRestoreLocalVector(rd->da,X0loc);CHKERRQ(ierr); 301 ierr = DMRestoreLocalVector(rd->da,Xloc);CHKERRQ(ierr); 302 ierr = DMRestoreLocalVector(rd->da,Xloc_t);CHKERRQ(ierr); 303 PetscFunctionReturn(0); 304 } 305 306 static PetscErrorCode RDCheckDomain_Private(RD rd,TS ts,Vec X,PetscBool *in) 307 { 308 PetscErrorCode ierr; 309 PetscInt minloc; 310 PetscReal min; 311 312 PetscFunctionBeginUser; 313 ierr = VecMin(X,&minloc,&min);CHKERRQ(ierr); 314 if (min < 0) { 315 SNES snes; 316 *in = PETSC_FALSE; 317 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 318 ierr = SNESSetFunctionDomainError(snes);CHKERRQ(ierr); 319 ierr = PetscInfo3(ts,"Domain violation at %D field %D value %g\n",minloc/2,minloc%2,(double)min);CHKERRQ(ierr); 320 } else *in = PETSC_TRUE; 321 PetscFunctionReturn(0); 322 } 323 324 /* Energy and temperature must remain positive */ 325 #define RDCheckDomain(rd,ts,X) do { \ 326 PetscErrorCode _ierr; \ 327 PetscBool _in; \ 328 _ierr = RDCheckDomain_Private(rd,ts,X,&_in);CHKERRQ(_ierr); \ 329 if (!_in) PetscFunctionReturn(0); \ 330 } while (0) 331 332 static PetscErrorCode RDIFunction_FD(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx) 333 { 334 PetscErrorCode ierr; 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; 341 342 PetscFunctionBeginUser; 343 RDCheckDomain(rd,ts,X); 344 ierr = RDGetLocalArrays(rd,ts,X,Xdot,&Theta,&dt,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);CHKERRQ(ierr); 345 ierr = DMDAVecGetArray(rd->da,F,&f);CHKERRQ(ierr); 346 ierr = DMDAGetLocalInfo(rd->da,&info);CHKERRQ(ierr); 347 VecZeroEntries(F); 348 349 hx = rd->L / (info.mx-1); 350 351 for (i=info.xs; i<info.xs+info.xm; i++) { 352 PetscReal rho = rd->rho; 353 PetscScalar Em_t,rad; 354 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); 369 370 if (i == 0) { /* Left boundary condition */ 371 PetscScalar D_R,bcTheta = rd->bcmidpoint ? Theta : 1.; 372 RDNode n, nx; 373 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 ierr = RDRestoreLocalArrays(rd,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);CHKERRQ(ierr); 396 ierr = DMDAVecRestoreArray(rd->da,F,&f);CHKERRQ(ierr); 397 if (rd->monitor_residual) {ierr = RDStateView(rd,X,Xdot,F);CHKERRQ(ierr);} 398 PetscFunctionReturn(0); 399 } 400 401 static PetscErrorCode RDIJacobian_FD(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) 402 { 403 PetscErrorCode ierr; 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; 410 411 PetscFunctionBeginUser; 412 RDCheckDomain(rd,ts,X); 413 ierr = RDGetLocalArrays(rd,ts,X,Xdot,&Theta,&dt,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);CHKERRQ(ierr); 414 ierr = DMDAGetLocalInfo(rd->da,&info);CHKERRQ(ierr); 415 hx = rd->L / (info.mx-1); 416 ierr = MatZeroEntries(B);CHKERRQ(ierr); 417 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; 423 424 /*rad = (1.-Theta)* */ RDRadiation(rd,&x0[i],0); /* + Theta* */ RDRadiation(rd,&x[i],&drad); 425 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; 439 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 } 451 452 ierr = PetscMemzero(K,sizeof(K));CHKERRQ(ierr); 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; 458 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 } 490 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); 493 494 col[0] = i-1; 495 col[1] = i; 496 col[2] = i+1<info.mx ? i+1 : -1; 497 ierr = MatSetValuesBlocked(B,1,&i,3,col,&K[0][0],INSERT_VALUES);CHKERRQ(ierr); 498 } 499 ierr = RDRestoreLocalArrays(rd,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);CHKERRQ(ierr); 500 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 501 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 502 if (A != B) { 503 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 504 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 505 } 506 PetscFunctionReturn(0); 507 } 508 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 } 521 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]; 529 530 PetscFunctionBeginUser; 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}; 551 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 } 564 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 PetscFunctionReturn(0); 573 } 574 575 /* 576 Finite element version 577 */ 578 static PetscErrorCode RDIFunction_FE(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx) 579 { 580 PetscErrorCode ierr; 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; 587 588 PetscFunctionBeginUser; 589 RDCheckDomain(rd,ts,X); 590 ierr = RDGetLocalArrays(rd,ts,X,Xdot,&Theta,&dt,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);CHKERRQ(ierr); 591 592 ierr = DMGetLocalVector(rd->da,&Floc);CHKERRQ(ierr); 593 ierr = VecZeroEntries(Floc);CHKERRQ(ierr); 594 ierr = DMDAVecGetArray(rd->da,Floc,&f);CHKERRQ(ierr); 595 ierr = DMDAGetLocalInfo(rd->da,&info);CHKERRQ(ierr); 596 597 /* Set up shape functions and quadrature for elements (assumes a uniform grid) */ 598 hx = rd->L / (info.mx-1); 599 ierr = RDGetQuadrature(rd,hx,&nq,weight,interp,deriv);CHKERRQ(ierr); 600 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); 609 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; 636 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 } 654 655 ierr = RDRestoreLocalArrays(rd,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);CHKERRQ(ierr); 656 ierr = DMDAVecRestoreArray(rd->da,Floc,&f);CHKERRQ(ierr); 657 ierr = VecZeroEntries(F);CHKERRQ(ierr); 658 ierr = DMLocalToGlobalBegin(rd->da,Floc,ADD_VALUES,F);CHKERRQ(ierr); 659 ierr = DMLocalToGlobalEnd(rd->da,Floc,ADD_VALUES,F);CHKERRQ(ierr); 660 ierr = DMRestoreLocalVector(rd->da,&Floc);CHKERRQ(ierr); 661 662 if (rd->monitor_residual) {ierr = RDStateView(rd,X,Xdot,F);CHKERRQ(ierr);} 663 PetscFunctionReturn(0); 664 } 665 666 static PetscErrorCode RDIJacobian_FE(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) 667 { 668 PetscErrorCode ierr; 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]; 676 677 PetscFunctionBeginUser; 678 RDCheckDomain(rd,ts,X); 679 ierr = RDGetLocalArrays(rd,ts,X,Xdot,&Theta,&dt,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);CHKERRQ(ierr); 680 ierr = DMDAGetLocalInfo(rd->da,&info);CHKERRQ(ierr); 681 hx = rd->L / (info.mx-1); 682 ierr = RDGetQuadrature(rd,hx,&nq,weight,interp,deriv);CHKERRQ(ierr); 683 ierr = MatZeroEntries(B);CHKERRQ(ierr); 684 for (i=info.xs; i<PetscMin(info.xs+info.xm,info.mx-1); i++) { 685 PetscInt rc[2]; 686 687 rc[0] = i; rc[1] = i+1; 688 ierr = PetscMemzero(K,sizeof(K));CHKERRQ(ierr); 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 ierr = MatSetValuesBlocked(B,2,rc,2,rc,&K[0][0],ADD_VALUES);CHKERRQ(ierr); 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; 717 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 ierr = MatSetValue(B,0,0,ratio*0.5,ADD_VALUES);CHKERRQ(ierr); 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 } 733 734 ierr = RDRestoreLocalArrays(rd,&X0loc,&x0,&Xloc,&x,&Xloc_t,&xdot);CHKERRQ(ierr); 735 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 736 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 737 if (A != B) { 738 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 739 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 740 } 741 PetscFunctionReturn(0); 742 } 743 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); } 746 747 static PetscErrorCode RDInitialState(RD rd,Vec X) 748 { 749 DMDALocalInfo info; 750 PetscInt i; 751 RDNode *x; 752 PetscErrorCode ierr; 753 754 PetscFunctionBeginUser; 755 ierr = DMDAGetLocalInfo(rd->da,&info);CHKERRQ(ierr); 756 ierr = DMDAVecGetArray(rd->da,X,&x);CHKERRQ(ierr); 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 ierr = DMDAVecRestoreArray(rd->da,X,&x);CHKERRQ(ierr); 776 PetscFunctionReturn(0); 777 } 778 779 static PetscErrorCode RDView(RD rd,Vec X,PetscViewer viewer) 780 { 781 PetscErrorCode ierr; 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; 789 790 PetscFunctionBeginUser; 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 ierr = DMDAGetInfo(rd->da,0, &M,0,0, 0,0,0, 0,0,0,0,0,0);CHKERRQ(ierr); 797 ierr = DMDAGetOwnershipRanges(rd->da,&lx,0,0);CHKERRQ(ierr); 798 ierr = PetscObjectGetComm((PetscObject)rd->da,&comm);CHKERRQ(ierr); 799 ierr = DMDACreate1d(comm,DM_BOUNDARY_NONE,M,1,0,lx,&da);CHKERRQ(ierr); 800 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 801 ierr = DMSetUp(da);CHKERRQ(ierr); 802 ierr = DMDASetUniformCoordinates(da,0.,rd->L,0.,0.,0.,0.);CHKERRQ(ierr); 803 ierr = DMDASetFieldName(da,0,"T_rad");CHKERRQ(ierr); 804 ierr = DMCreateGlobalVector(da,&Y);CHKERRQ(ierr); 805 806 /* Compute the radiation temperature from the solution at each node */ 807 ierr = VecGetLocalSize(Y,&m);CHKERRQ(ierr); 808 ierr = VecGetArrayRead(X,(const PetscScalar **)&x);CHKERRQ(ierr); 809 ierr = VecGetArray(Y,&y);CHKERRQ(ierr); 810 for (i=0; i<m; i++) y[i] = RDRadiationTemperature(rd,x[i].E); 811 ierr = VecRestoreArrayRead(X,(const PetscScalar**)&x);CHKERRQ(ierr); 812 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr); 813 814 ierr = VecView(Y,viewer);CHKERRQ(ierr); 815 ierr = VecDestroy(&Y);CHKERRQ(ierr); 816 ierr = DMDestroy(&da);CHKERRQ(ierr); 817 PetscFunctionReturn(0); 818 } 819 820 static PetscErrorCode RDTestDifferentiation(RD rd) 821 { 822 MPI_Comm comm; 823 PetscErrorCode ierr; 824 RDNode n,nx; 825 PetscScalar epsilon; 826 827 PetscFunctionBeginUser; 828 ierr = PetscObjectGetComm((PetscObject)rd->da,&comm);CHKERRQ(ierr); 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 ierr = 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));CHKERRQ(ierr); 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 ierr = 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));CHKERRQ(ierr); 862 ierr = 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));CHKERRQ(ierr); 864 } 865 { 866 PetscInt i; 867 PetscReal hx = 1.; 868 PetscScalar a0; 869 RDNode n0[3],n1[3],d[3],fd[3]; 870 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 ierr = PetscMemcpy(n1,n0,sizeof(n0));CHKERRQ(ierr); n1[i].E += epsilon; 880 fd[i].E = (RDDiffusion(rd,hx,n1,1,0)-a0)/epsilon; 881 ierr = PetscMemcpy(n1,n0,sizeof(n0));CHKERRQ(ierr); n1[i].T += epsilon; 882 fd[i].T = (RDDiffusion(rd,hx,n1,1,0)-a0)/epsilon; 883 ierr = 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));CHKERRQ(ierr); 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 ierr = 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));CHKERRQ(ierr); 898 } 899 PetscFunctionReturn(0); 900 } 901 902 static PetscErrorCode RDCreate(MPI_Comm comm,RD *inrd) 903 { 904 PetscErrorCode ierr; 905 RD rd; 906 PetscReal meter=0,kilogram=0,second=0,Kelvin=0,Joule=0,Watt=0; 907 908 PetscFunctionBeginUser; 909 *inrd = 0; 910 ierr = PetscNew(&rd);CHKERRQ(ierr); 911 912 ierr = PetscOptionsBegin(comm,NULL,"Options for nonequilibrium radiation-diffusion with RD ionization",NULL);CHKERRQ(ierr); 913 { 914 rd->initial = 1; 915 ierr = PetscOptionsInt("-rd_initial","Initial condition (1=Marshak, 2=Blast, 3=Marshak+)","",rd->initial,&rd->initial,0);CHKERRQ(ierr); 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 ierr = PetscOptionsReal("-rd_unit_meter","Length of 1 meter in nondimensional units","",rd->unit.meter,&rd->unit.meter,0);CHKERRQ(ierr); 934 ierr = PetscOptionsReal("-rd_unit_kilogram","Mass of 1 kilogram in nondimensional units","",rd->unit.kilogram,&rd->unit.kilogram,0);CHKERRQ(ierr); 935 ierr = PetscOptionsReal("-rd_unit_second","Time of a second in nondimensional units","",rd->unit.second,&rd->unit.second,0);CHKERRQ(ierr); 936 ierr = PetscOptionsReal("-rd_unit_Kelvin","Temperature of a Kelvin in nondimensional units","",rd->unit.Kelvin,&rd->unit.Kelvin,0);CHKERRQ(ierr); 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; 947 948 ierr = PetscOptionsBool("-rd_monitor_residual","Display residuals every time they are evaluated","",rd->monitor_residual,&rd->monitor_residual,NULL);CHKERRQ(ierr); 949 ierr = PetscOptionsEnum("-rd_discretization","Discretization type","",DiscretizationTypes,(PetscEnum)rd->discretization,(PetscEnum*)&rd->discretization,NULL);CHKERRQ(ierr); 950 if (rd->discretization == DISCRETIZATION_FE) { 951 rd->quadrature = QUADRATURE_GAUSS2; 952 ierr = PetscOptionsEnum("-rd_quadrature","Finite element quadrature","",QuadratureTypes,(PetscEnum)rd->quadrature,(PetscEnum*)&rd->quadrature,NULL);CHKERRQ(ierr); 953 } 954 ierr = PetscOptionsEnum("-rd_jacobian","Type of finite difference Jacobian","",JacobianTypes,(PetscEnum)rd->jacobian,(PetscEnum*)&rd->jacobian,NULL);CHKERRQ(ierr); 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 ierr = PetscOptionsEnum("-rd_leftbc","Left boundary condition","",BCTypes,(PetscEnum)rd->leftbc,(PetscEnum*)&rd->leftbc,NULL);CHKERRQ(ierr); 983 ierr = PetscOptionsReal("-rd_E_applied","Radiation flux at left end of domain","",rd->Eapplied,&rd->Eapplied,NULL);CHKERRQ(ierr); 984 ierr = PetscOptionsReal("-rd_beta","Thermal exponent for photon absorption","",rd->beta,&rd->beta,NULL);CHKERRQ(ierr); 985 ierr = PetscOptionsReal("-rd_gamma","Thermal exponent for diffusion coefficient","",rd->gamma,&rd->gamma,NULL);CHKERRQ(ierr); 986 ierr = PetscOptionsBool("-rd_view_draw","Draw final solution","",rd->view_draw,&rd->view_draw,NULL);CHKERRQ(ierr); 987 ierr = PetscOptionsBool("-rd_endpoint","Discretize using endpoints (like trapezoid rule) instead of midpoint","",rd->endpoint,&rd->endpoint,NULL);CHKERRQ(ierr); 988 ierr = PetscOptionsBool("-rd_bcmidpoint","Impose the boundary condition at the midpoint (Theta) of the interval","",rd->bcmidpoint,&rd->bcmidpoint,NULL);CHKERRQ(ierr); 989 ierr = PetscOptionsBool("-rd_bclimit","Limit diffusion coefficient in definition of Robin boundary condition","",rd->bclimit,&rd->bclimit,NULL);CHKERRQ(ierr); 990 ierr = PetscOptionsBool("-rd_test_diff","Test differentiation in constitutive relations","",rd->test_diff,&rd->test_diff,NULL);CHKERRQ(ierr); 991 ierr = PetscOptionsString("-rd_view_binary","File name to hold final solution","",rd->view_binary,rd->view_binary,sizeof(rd->view_binary),NULL);CHKERRQ(ierr); 992 } 993 ierr = PetscOptionsEnd();CHKERRQ(ierr); 994 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 } 1020 1021 ierr = DMDACreate1d(comm,DM_BOUNDARY_NONE,20,sizeof(RDNode)/sizeof(PetscScalar),1,NULL,&rd->da);CHKERRQ(ierr); 1022 ierr = DMSetFromOptions(rd->da);CHKERRQ(ierr); 1023 ierr = DMSetUp(rd->da);CHKERRQ(ierr); 1024 ierr = DMDASetFieldName(rd->da,0,"E");CHKERRQ(ierr); 1025 ierr = DMDASetFieldName(rd->da,1,"T");CHKERRQ(ierr); 1026 ierr = DMDASetUniformCoordinates(rd->da,0.,1.,0.,0.,0.,0.);CHKERRQ(ierr); 1027 1028 *inrd = rd; 1029 PetscFunctionReturn(0); 1030 } 1031 1032 int main(int argc, char *argv[]) 1033 { 1034 PetscErrorCode ierr; 1035 RD rd; 1036 TS ts; 1037 SNES snes; 1038 Vec X; 1039 Mat A,B; 1040 PetscInt steps; 1041 PetscReal ftime; 1042 1043 ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr; 1044 ierr = RDCreate(PETSC_COMM_WORLD,&rd);CHKERRQ(ierr); 1045 ierr = DMCreateGlobalVector(rd->da,&X);CHKERRQ(ierr); 1046 ierr = DMSetMatType(rd->da,MATAIJ);CHKERRQ(ierr); 1047 ierr = DMCreateMatrix(rd->da,&B);CHKERRQ(ierr); 1048 ierr = RDInitialState(rd,X);CHKERRQ(ierr); 1049 1050 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 1051 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 1052 ierr = TSSetType(ts,TSTHETA);CHKERRQ(ierr); 1053 ierr = TSSetDM(ts,rd->da);CHKERRQ(ierr); 1054 switch (rd->discretization) { 1055 case DISCRETIZATION_FD: 1056 ierr = TSSetIFunction(ts,NULL,RDIFunction_FD,rd);CHKERRQ(ierr); 1057 if (rd->jacobian == JACOBIAN_ANALYTIC) ierr = TSSetIJacobian(ts,B,B,RDIJacobian_FD,rd);CHKERRQ(ierr); 1058 break; 1059 case DISCRETIZATION_FE: 1060 ierr = TSSetIFunction(ts,NULL,RDIFunction_FE,rd);CHKERRQ(ierr); 1061 if (rd->jacobian == JACOBIAN_ANALYTIC) ierr = TSSetIJacobian(ts,B,B,RDIJacobian_FE,rd);CHKERRQ(ierr); 1062 break; 1063 } 1064 ierr = TSSetMaxTime(ts,rd->final_time);CHKERRQ(ierr); 1065 ierr = TSSetTimeStep(ts,1e-3);CHKERRQ(ierr); 1066 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 1067 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 1068 1069 A = B; 1070 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 1071 switch (rd->jacobian) { 1072 case JACOBIAN_ANALYTIC: 1073 break; 1074 case JACOBIAN_MATRIXFREE: 1075 break; 1076 case JACOBIAN_FD_COLORING: { 1077 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobianDefaultColor,0);CHKERRQ(ierr); 1078 } break; 1079 case JACOBIAN_FD_FULL: 1080 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobianDefault,ts);CHKERRQ(ierr); 1081 break; 1082 } 1083 1084 if (rd->test_diff) { 1085 ierr = RDTestDifferentiation(rd);CHKERRQ(ierr); 1086 } 1087 ierr = TSSolve(ts,X);CHKERRQ(ierr); 1088 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 1089 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 1090 ierr = PetscPrintf(PETSC_COMM_WORLD,"Steps %D final time %g\n",steps,(double)ftime);CHKERRQ(ierr); 1091 if (rd->view_draw) { 1092 ierr = RDView(rd,X,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 1093 } 1094 if (rd->view_binary[0]) { 1095 PetscViewer viewer; 1096 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,rd->view_binary,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 1097 ierr = RDView(rd,X,viewer);CHKERRQ(ierr); 1098 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 1099 } 1100 ierr = VecDestroy(&X);CHKERRQ(ierr); 1101 ierr = MatDestroy(&B);CHKERRQ(ierr); 1102 ierr = RDDestroy(&rd);CHKERRQ(ierr); 1103 ierr = TSDestroy(&ts);CHKERRQ(ierr); 1104 ierr = PetscFinalize(); 1105 return ierr; 1106 } 1107 /*TEST 1108 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 1112 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 1117 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 1123 1124 TEST*/ 1125