Lines Matching refs:rd
82 static PetscErrorCode RDDestroy(RD *rd) in RDDestroy() argument
85 PetscCall(DMDestroy(&(*rd)->da)); in RDDestroy()
86 PetscCall(PetscFree(*rd)); in RDDestroy()
98 static void RDMaterialEnergy(RD rd, const RDNode *n, PetscScalar *Em, RDNode *dEm) in RDMaterialEnergy() argument
100 rd->MaterialEnergy(rd, n, Em, dEm); in RDMaterialEnergy()
113 static void RDMaterialEnergy_Saha(RD rd, const RDNode *n, PetscScalar *inEm, RDNode *dEm) in RDMaterialEnergy_Saha() argument
115 …rd->I_H / (rd->k * T), chi_t = -chi / T * T_t, a = 1., a_t = 0, b = 4. * rd->m_p / rd->rho * Petsc… in RDMaterialEnergy_Saha()
118 Em = rd->k * T / rd->m_p * (1.5 * (1. + alpha) + alpha * chi); /* Eq 6 */ in RDMaterialEnergy_Saha()
122 dEm->T = Em / T * T_t + rd->k * T / rd->m_p * (1.5 * alpha_t + alpha_t * chi + alpha * chi_t); in RDMaterialEnergy_Saha()
126 static void RDMaterialEnergy_Reduced(RD rd, const RDNode *n, PetscScalar *Em, RDNode *dEm) in RDMaterialEnergy_Reduced() argument
138 static void RDSigma_R(RD rd, RDNode *n, PetscScalar *sigma_R, RDNode *dsigma_R) in RDSigma_R() argument
140 *sigma_R = rd->K_R * rd->rho * PetscPowScalar(n->T, -rd->gamma); in RDSigma_R()
142 dsigma_R->T = -rd->gamma * (*sigma_R) / n->T; in RDSigma_R()
146 static void RDDiffusionCoefficient(RD rd, PetscBool limit, RDNode *n, RDNode *nx, PetscScalar *D_R,… in RDDiffusionCoefficient() argument
151 RDSigma_R(rd, n, &sigma_R, &dsigma_R); in RDDiffusionCoefficient()
152 denom = 3. * rd->rho * sigma_R + (int)limit * PetscAbsScalar(nx->E) / n->E; in RDDiffusionCoefficient()
154 ddenom.T = 3. * rd->rho * dsigma_R.T; in RDDiffusionCoefficient()
157 *D_R = rd->c / denom; in RDDiffusionCoefficient()
159 dD_R->E = -rd->c / PetscSqr(denom) * ddenom.E; in RDDiffusionCoefficient()
160 dD_R->T = -rd->c / PetscSqr(denom) * ddenom.T; in RDDiffusionCoefficient()
163 dxD_R->E = -rd->c / PetscSqr(denom) * dxdenom.E; in RDDiffusionCoefficient()
164 dxD_R->T = -rd->c / PetscSqr(denom) * dxdenom.T; in RDDiffusionCoefficient()
168 static PetscErrorCode RDStateView(RD rd, Vec X, Vec Xdot, Vec F) in RDStateView() argument
176 PetscCall(PetscObjectGetComm((PetscObject)rd->da, &comm)); in RDStateView()
177 PetscCall(DMDAGetLocalInfo(rd->da, &info)); in RDStateView()
178 PetscCall(DMDAVecGetArrayRead(rd->da, X, (void *)&x)); in RDStateView()
179 PetscCall(DMDAVecGetArrayRead(rd->da, Xdot, (void *)&xdot)); in RDStateView()
180 PetscCall(DMDAVecGetArrayRead(rd->da, F, (void *)&f)); in RDStateView()
185 PetscCall(DMDAVecRestoreArrayRead(rd->da, X, (void *)&x)); in RDStateView()
186 PetscCall(DMDAVecRestoreArrayRead(rd->da, Xdot, (void *)&xdot)); in RDStateView()
187 PetscCall(DMDAVecRestoreArrayRead(rd->da, F, (void *)&f)); in RDStateView()
192 static PetscScalar RDRadiation(RD rd, const RDNode *n, RDNode *dn) in RDRadiation() argument
194 …rd->K_p * rd->rho * PetscPowScalar(n->T, -rd->beta), sigma_p_T = -rd->beta * sigma_p / n->T, tmp =… in RDRadiation()
202 static PetscScalar RDDiffusion(RD rd, PetscReal hx, const RDNode x[], PetscInt i, RDNode d[]) in RDDiffusion() argument
212 RDDiffusionCoefficient(rd, PETSC_TRUE, &n_L, &nx_L, &D_L, &dD_L, &dxD_L); in RDDiffusion()
223 RDDiffusionCoefficient(rd, PETSC_TRUE, &n_R, &nx_R, &D_R, &dD_R, &dxD_R); in RDDiffusion()
241 static PetscErrorCode RDGetLocalArrays(RD rd, TS ts, Vec X, Vec Xdot, PetscReal *Theta, PetscReal *… in RDGetLocalArrays() argument
246 PetscCall(DMGetLocalVector(rd->da, X0loc)); in RDGetLocalArrays()
247 PetscCall(DMGetLocalVector(rd->da, Xloc)); in RDGetLocalArrays()
248 PetscCall(DMGetLocalVector(rd->da, Xloc_t)); in RDGetLocalArrays()
250 PetscCall(DMGlobalToLocalBegin(rd->da, X, INSERT_VALUES, *Xloc)); in RDGetLocalArrays()
251 PetscCall(DMGlobalToLocalEnd(rd->da, X, INSERT_VALUES, *Xloc)); in RDGetLocalArrays()
252 PetscCall(DMGlobalToLocalBegin(rd->da, Xdot, INSERT_VALUES, *Xloc_t)); in RDGetLocalArrays()
253 PetscCall(DMGlobalToLocalEnd(rd->da, Xdot, INSERT_VALUES, *Xloc_t)); in RDGetLocalArrays()
261 if (istheta && rd->endpoint) PetscCall(TSThetaGetTheta(ts, Theta)); in RDGetLocalArrays()
266 …if (rd->endpoint) PetscCall(VecWAXPY(*Xloc, *dt, *Xloc_t, *X0loc)); /* move the abscissa to the en… in RDGetLocalArrays()
268 PetscCall(DMDAVecGetArray(rd->da, *X0loc, x0)); in RDGetLocalArrays()
269 PetscCall(DMDAVecGetArray(rd->da, *Xloc, x)); in RDGetLocalArrays()
270 PetscCall(DMDAVecGetArray(rd->da, *Xloc_t, xdot)); in RDGetLocalArrays()
274 static PetscErrorCode RDRestoreLocalArrays(RD rd, Vec *X0loc, RDNode **x0, Vec *Xloc, RDNode **x, V… in RDRestoreLocalArrays() argument
277 PetscCall(DMDAVecRestoreArray(rd->da, *X0loc, x0)); in RDRestoreLocalArrays()
278 PetscCall(DMDAVecRestoreArray(rd->da, *Xloc, x)); in RDRestoreLocalArrays()
279 PetscCall(DMDAVecRestoreArray(rd->da, *Xloc_t, xdot)); in RDRestoreLocalArrays()
280 PetscCall(DMRestoreLocalVector(rd->da, X0loc)); in RDRestoreLocalArrays()
281 PetscCall(DMRestoreLocalVector(rd->da, Xloc)); in RDRestoreLocalArrays()
282 PetscCall(DMRestoreLocalVector(rd->da, Xloc_t)); in RDRestoreLocalArrays()
286 static PetscErrorCode PETSC_UNUSED RDCheckDomain_Private(RD rd, TS ts, Vec X, PetscBool *in) in RDCheckDomain_Private() argument
304 #define RDCheckDomain(rd, ts, X) \ argument
307 PetscCall(RDCheckDomain_Private(rd, ts, X, &_in)); \
313 RD rd = (RD)ctx; in RDIFunction_FD() local
321 PetscCall(RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); in RDIFunction_FD()
322 PetscCall(DMDAVecGetArray(rd->da, F, &f)); in RDIFunction_FD()
323 PetscCall(DMDAGetLocalInfo(rd->da, &info)); in RDIFunction_FD()
326 hx = rd->L / (info.mx - 1); in RDIFunction_FD()
329 PetscReal rho = rd->rho; in RDIFunction_FD()
332 rad = (1. - Theta) * RDRadiation(rd, &x0[i], 0) + Theta * RDRadiation(rd, &x[i], 0); in RDIFunction_FD()
333 if (rd->endpoint) { in RDIFunction_FD()
335 RDMaterialEnergy(rd, &x0[i], &Em0, NULL); in RDIFunction_FD()
336 RDMaterialEnergy(rd, &x[i], &Em1, NULL); in RDIFunction_FD()
340 RDMaterialEnergy(rd, &x[i], NULL, &dEm); in RDIFunction_FD()
348 PetscScalar D_R, bcTheta = rd->bcmidpoint ? Theta : 1.; in RDIFunction_FD()
355 switch (rd->leftbc) { in RDIFunction_FD()
357 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D_R, 0, 0); in RDIFunction_FD()
358 f[0].E = hx * (n.E - 2. * D_R * nx.E - rd->Eapplied); in RDIFunction_FD()
364 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT, rd->initial); in RDIFunction_FD()
369 …PetscScalar diff = (1. - Theta) * RDDiffusion(rd, hx, x0, i, 0) + Theta * RDDiffusion(rd, hx, x, i… in RDIFunction_FD()
373 PetscCall(RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); in RDIFunction_FD()
374 PetscCall(DMDAVecRestoreArray(rd->da, F, &f)); in RDIFunction_FD()
375 if (rd->monitor_residual) PetscCall(RDStateView(rd, X, Xdot, F)); in RDIFunction_FD()
381 RD rd = (RD)ctx; in RDIJacobian_FD() local
389 PetscCall(RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); in RDIJacobian_FD()
390 PetscCall(DMDAGetLocalInfo(rd->da, &info)); in RDIJacobian_FD()
391 hx = rd->L / (info.mx - 1); in RDIJacobian_FD()
396 PetscReal rho = rd->rho; in RDIJacobian_FD()
400 /*rad = (1.-Theta)* */ RDRadiation(rd, &x0[i], 0); /* + Theta* */ in RDIJacobian_FD()
401 RDRadiation(rd, &x[i], &drad); in RDIJacobian_FD()
403 if (rd->endpoint) { in RDIJacobian_FD()
406 RDMaterialEnergy(rd, &x0[i], &Em0, NULL); in RDIJacobian_FD()
407 RDMaterialEnergy(rd, &x[i], &Em1, &dEm1); in RDIJacobian_FD()
419 RDMaterialEnergy(rd, &x[i], NULL, &dEm); in RDIJacobian_FD()
420 RDMaterialEnergy(rd, &n1, NULL, &dEm1); in RDIJacobian_FD()
432 PetscScalar D, bcTheta = rd->bcmidpoint ? Theta : 1.; in RDIJacobian_FD()
440 switch (rd->leftbc) { in RDIJacobian_FD()
442 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, &dD, &dxD); in RDIJacobian_FD()
452 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT, rd->initial); in RDIJacobian_FD()
460 /*diff = (1.-Theta)*RDDiffusion(rd,hx,x0,i,0) + Theta* */ RDDiffusion(rd, hx, x, i, ddiff); in RDIJacobian_FD()
477 PetscCall(RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); in RDIJacobian_FD()
506 static PetscErrorCode RDGetQuadrature(RD rd, PetscReal hx, PetscInt *nq, PetscReal weight[], PetscR… in RDGetQuadrature() argument
512 switch (rd->quadrature) { in RDGetQuadrature()
611 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown quadrature %d", (int)rd->quadrature); in RDGetQuadrature()
629 RD rd = (RD)ctx; in RDIFunction_FE() local
637 PetscCall(RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); in RDIFunction_FE()
639 PetscCall(DMGetLocalVector(rd->da, &Floc)); in RDIFunction_FE()
641 PetscCall(DMDAVecGetArray(rd->da, Floc, &f)); in RDIFunction_FE()
642 PetscCall(DMDAGetLocalInfo(rd->da, &info)); in RDIFunction_FE()
645 hx = rd->L / (info.mx - 1); in RDIFunction_FE()
646 PetscCall(RDGetQuadrature(rd, hx, &nq, weight, interp, deriv)); in RDIFunction_FE()
650 PetscReal rho = rd->rho; in RDIFunction_FE()
657 rad = (1. - Theta) * RDRadiation(rd, &n0, 0) + Theta * RDRadiation(rd, &n, 0); in RDIFunction_FE()
658 if (rd->endpoint) { in RDIFunction_FE()
660 RDMaterialEnergy(rd, &n0, &Em0, NULL); in RDIFunction_FE()
661 RDMaterialEnergy(rd, &n, &Em1, NULL); in RDIFunction_FE()
665 RDMaterialEnergy(rd, &n, NULL, &dEm); in RDIFunction_FE()
668 RDDiffusionCoefficient(rd, PETSC_TRUE, &n0, &n0x, &D0_R, 0, 0); in RDIFunction_FE()
669 RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, 0, 0); in RDIFunction_FE()
677 switch (rd->leftbc) { in RDIFunction_FE()
680 PetscReal ratio, bcTheta = rd->bcmidpoint ? Theta : 1.; in RDIFunction_FE()
687 RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, 0, 0); in RDIFunction_FE()
688 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D_R_bc, 0, 0); in RDIFunction_FE()
692 f[0].E += -ratio * 0.5 * (rd->Eapplied - n.E); in RDIFunction_FE()
698 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT, rd->initial); in RDIFunction_FE()
702 PetscCall(RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); in RDIFunction_FE()
703 PetscCall(DMDAVecRestoreArray(rd->da, Floc, &f)); in RDIFunction_FE()
705 PetscCall(DMLocalToGlobalBegin(rd->da, Floc, ADD_VALUES, F)); in RDIFunction_FE()
706 PetscCall(DMLocalToGlobalEnd(rd->da, Floc, ADD_VALUES, F)); in RDIFunction_FE()
707 PetscCall(DMRestoreLocalVector(rd->da, &Floc)); in RDIFunction_FE()
709 if (rd->monitor_residual) PetscCall(RDStateView(rd, X, Xdot, F)); in RDIFunction_FE()
715 RD rd = (RD)ctx; in RDIJacobian_FE() local
724 PetscCall(RDGetLocalArrays(rd, ts, X, Xdot, &Theta, &dt, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); in RDIJacobian_FE()
725 PetscCall(DMDAGetLocalInfo(rd->da, &info)); in RDIJacobian_FE()
726 hx = rd->L / (info.mx - 1); in RDIJacobian_FE()
727 PetscCall(RDGetQuadrature(rd, hx, &nq, weight, interp, deriv)); in RDIJacobian_FE()
741 rad = RDRadiation(rd, &n, &drad); in RDIJacobian_FE()
742 RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, &dD_R, &dxD_R); in RDIJacobian_FE()
743 RDMaterialEnergy(rd, &n, NULL, &dEm); in RDIJacobian_FE()
749 …K[j * 2 + 1][k * 2 + 1] += interp[q][j] * weight[q] * (a * rd->rho * dEm.T + drad.T) * interp[q][k… in RDIJacobian_FE()
756 switch (rd->leftbc) { in RDIJacobian_FE()
766 RDDiffusionCoefficient(rd, PETSC_TRUE, &n, &nx, &D_R, 0, 0); in RDIJacobian_FE()
767 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D_R_bc, 0, 0); in RDIJacobian_FE()
775 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT, rd->initial); in RDIJacobian_FE()
779 PetscCall(RDRestoreLocalArrays(rd, &X0loc, &x0, &Xloc, &x, &Xloc_t, &xdot)); in RDIJacobian_FE()
790 static PetscScalar RDRadiationTemperature(RD rd, PetscScalar E) in RDRadiationTemperature() argument
792 return PetscPowScalar(E * rd->c / (4. * rd->sigma_b), 0.25); in RDRadiationTemperature()
795 static PetscErrorCode RDInitialState(RD rd, Vec X) in RDInitialState() argument
802 PetscCall(DMDAGetLocalInfo(rd->da, &info)); in RDInitialState()
803 PetscCall(DMDAVecGetArray(rd->da, X, &x)); in RDInitialState()
805 PetscReal coord = i * rd->L / (info.mx - 1); in RDInitialState()
806 switch (rd->initial) { in RDInitialState()
809 x[i].T = RDRadiationTemperature(rd, x[i].E); in RDInitialState()
813 x[i].T = RDRadiationTemperature(rd, x[i].E); in RDInitialState()
816 x[i].E = 7.56e-2 * rd->unit.Joule / PetscPowScalarInt(rd->unit.meter, 3); in RDInitialState()
817 x[i].T = RDRadiationTemperature(rd, x[i].E); in RDInitialState()
820 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No initial state %" PetscInt_FMT, rd->initial); in RDInitialState()
823 PetscCall(DMDAVecRestoreArray(rd->da, X, &x)); in RDInitialState()
827 static PetscErrorCode RDView(RD rd, Vec X, PetscViewer viewer) in RDView() argument
843 PetscCall(DMDAGetInfo(rd->da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); in RDView()
844 PetscCall(DMDAGetOwnershipRanges(rd->da, &lx, 0, 0)); in RDView()
845 PetscCall(PetscObjectGetComm((PetscObject)rd->da, &comm)); in RDView()
849 PetscCall(DMDASetUniformCoordinates(da, 0., rd->L, 0., 0., 0., 0.)); in RDView()
857 for (i = 0; i < m; i++) y[i] = RDRadiationTemperature(rd, x[i].E); in RDView()
867 static PetscErrorCode RDTestDifferentiation(RD rd) in RDTestDifferentiation() argument
874 PetscCall(PetscObjectGetComm((PetscObject)rd->da, &comm)); in RDTestDifferentiation()
881 rd->MaterialEnergy(rd, &n, &Em0, &dEm); in RDTestDifferentiation()
884 rd->MaterialEnergy(rd, &n, &Em1, 0); in RDTestDifferentiation()
888 rd->MaterialEnergy(rd, &n, &Em1, 0); in RDTestDifferentiation()
900 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D0, &dD, &dxD); in RDTestDifferentiation()
905 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0); in RDTestDifferentiation()
911 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0); in RDTestDifferentiation()
917 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0); in RDTestDifferentiation()
923 RDDiffusionCoefficient(rd, rd->bclimit, &n, &nx, &D, 0, 0); in RDTestDifferentiation()
942 a0 = RDDiffusion(rd, hx, n0, 1, d); in RDTestDifferentiation()
946 fd[i].E = (RDDiffusion(rd, hx, n1, 1, 0) - a0) / epsilon; in RDTestDifferentiation()
949 fd[i].T = (RDDiffusion(rd, hx, n1, 1, 0) - a0) / epsilon; in RDTestDifferentiation()
959 rad0 = RDRadiation(rd, &n, &drad); in RDTestDifferentiation()
962 rad = RDRadiation(rd, &n, 0); in RDTestDifferentiation()
966 rad = RDRadiation(rd, &n, 0); in RDTestDifferentiation()
976 RD rd; in RDCreate() local
981 PetscCall(PetscNew(&rd)); in RDCreate()
985 rd->initial = 1; in RDCreate()
986 …initial", "Initial condition (1=Marshak, 2=Blast, 3=Marshak+)", "", rd->initial, &rd->initial, 0)); in RDCreate()
987 switch (rd->initial) { in RDCreate()
990 rd->unit.kilogram = 1.; in RDCreate()
991 rd->unit.meter = 1.; in RDCreate()
992 rd->unit.second = 1.; in RDCreate()
993 rd->unit.Kelvin = 1.; in RDCreate()
996 rd->unit.kilogram = 1.e12; in RDCreate()
997 rd->unit.meter = 1.; in RDCreate()
998 rd->unit.second = 1.e9; in RDCreate()
999 rd->unit.Kelvin = 1.; in RDCreate()
1002 … SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown initial condition %" PetscInt_FMT, rd->initial); in RDCreate()
1005 …"-rd_unit_meter", "Length of 1 meter in nondimensional units", "", rd->unit.meter, &rd->unit.meter… in RDCreate()
1006 …nit_kilogram", "Mass of 1 kilogram in nondimensional units", "", rd->unit.kilogram, &rd->unit.kilo… in RDCreate()
1007 …-rd_unit_second", "Time of a second in nondimensional units", "", rd->unit.second, &rd->unit.secon… in RDCreate()
1008 …t_Kelvin", "Temperature of a Kelvin in nondimensional units", "", rd->unit.Kelvin, &rd->unit.Kelvi… in RDCreate()
1010 rd->unit.Joule = rd->unit.kilogram * PetscSqr(rd->unit.meter / rd->unit.second); in RDCreate()
1011 rd->unit.Watt = rd->unit.Joule / rd->unit.second; in RDCreate()
1013 meter = rd->unit.meter; in RDCreate()
1014 kilogram = rd->unit.kilogram; in RDCreate()
1015 second = rd->unit.second; in RDCreate()
1016 Kelvin = rd->unit.Kelvin; in RDCreate()
1017 Joule = rd->unit.Joule; in RDCreate()
1018 Watt = rd->unit.Watt; in RDCreate()
1020 …sidual", "Display residuals every time they are evaluated", "", rd->monitor_residual, &rd->monitor… in RDCreate()
1021 …Discretization type", "", DiscretizationTypes, (PetscEnum)rd->discretization, (PetscEnum *)&rd->di… in RDCreate()
1022 if (rd->discretization == DISCRETIZATION_FE) { in RDCreate()
1023 rd->quadrature = QUADRATURE_GAUSS2; in RDCreate()
1024 …Finite element quadrature", "", QuadratureTypes, (PetscEnum)rd->quadrature, (PetscEnum *)&rd->quad… in RDCreate()
1026 …f finite difference Jacobian", "", JacobianTypes, (PetscEnum)rd->jacobian, (PetscEnum *)&rd->jacob… in RDCreate()
1027 switch (rd->initial) { in RDCreate()
1029 rd->leftbc = BC_ROBIN; in RDCreate()
1030 rd->Eapplied = 4 * rd->unit.Joule / PetscPowRealInt(rd->unit.meter, 3); in RDCreate()
1031 rd->L = 1. * rd->unit.meter; in RDCreate()
1032 rd->beta = 3.0; in RDCreate()
1033 rd->gamma = 3.0; in RDCreate()
1034 rd->final_time = 3 * second; in RDCreate()
1037 rd->leftbc = BC_NEUMANN; in RDCreate()
1038 rd->Eapplied = 0.; in RDCreate()
1039 rd->L = 1. * rd->unit.meter; in RDCreate()
1040 rd->beta = 3.0; in RDCreate()
1041 rd->gamma = 3.0; in RDCreate()
1042 rd->final_time = 1 * second; in RDCreate()
1045 rd->leftbc = BC_ROBIN; in RDCreate()
1046 rd->Eapplied = 7.503e6 * rd->unit.Joule / PetscPowRealInt(rd->unit.meter, 3); in RDCreate()
1047 rd->L = 5. * rd->unit.meter; in RDCreate()
1048 rd->beta = 3.5; in RDCreate()
1049 rd->gamma = 3.5; in RDCreate()
1050 rd->final_time = 20e-9 * second; in RDCreate()
1053 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Initial %" PetscInt_FMT, rd->initial); in RDCreate()
1055 …d_leftbc", "Left boundary condition", "", BCTypes, (PetscEnum)rd->leftbc, (PetscEnum *)&rd->leftbc… in RDCreate()
1056 …nsReal("-rd_E_applied", "Radiation flux at left end of domain", "", rd->Eapplied, &rd->Eapplied, N… in RDCreate()
1057 …cOptionsReal("-rd_beta", "Thermal exponent for photon absorption", "", rd->beta, &rd->beta, NULL)); in RDCreate()
1058 …sReal("-rd_gamma", "Thermal exponent for diffusion coefficient", "", rd->gamma, &rd->gamma, NULL)); in RDCreate()
1059 …PetscCall(PetscOptionsBool("-rd_view_draw", "Draw final solution", "", rd->view_draw, &rd->view_dr… in RDCreate()
1060 …ize using endpoints (like trapezoid rule) instead of midpoint", "", rd->endpoint, &rd->endpoint, N… in RDCreate()
1061 …e boundary condition at the midpoint (Theta) of the interval", "", rd->bcmidpoint, &rd->bcmidpoint… in RDCreate()
1062 …ffusion coefficient in definition of Robin boundary condition", "", rd->bclimit, &rd->bclimit, NUL… in RDCreate()
1063 …_test_diff", "Test differentiation in constitutive relations", "", rd->test_diff, &rd->test_diff, … in RDCreate()
1064 …_view_binary", "File name to hold final solution", "", rd->view_binary, rd->view_binary, sizeof(rd… in RDCreate()
1068 switch (rd->initial) { in RDCreate()
1071 rd->rho = 1.; in RDCreate()
1072 rd->c = 1.; in RDCreate()
1073 rd->K_R = 1.; in RDCreate()
1074 rd->K_p = 1.; in RDCreate()
1075 rd->sigma_b = 0.25; in RDCreate()
1076 rd->MaterialEnergy = RDMaterialEnergy_Reduced; in RDCreate()
1080 rd->rho = 1.17e-3 * kilogram / (meter * meter * meter); /* density */ in RDCreate()
1081 …rd->K_R = 7.44e18 * PetscPowRealInt(meter, 5) * PetscPowReal(Kelvin, 3.5) * PetscPowRea… in RDCreate()
1082 …rd->K_p = 2.33e20 * PetscPowRealInt(meter, 5) * PetscPowReal(Kelvin, 3.5) * PetscPowRea… in RDCreate()
1083 …rd->I_H = 2.179e-18 * Joule; /*… in RDCreate()
1084 …rd->m_p = 1.673e-27 * kilogram; /*… in RDCreate()
1085 …rd->m_e = 9.109e-31 * kilogram; /*… in RDCreate()
1086 …rd->h = 6.626e-34 * Joule * second; /*… in RDCreate()
1087 …rd->k = 1.381e-23 * Joule / Kelvin; /*… in RDCreate()
1088 …rd->c = 3.00e8 * meter / second; /*… in RDCreate()
1089 …rd->sigma_b = 5.67e-8 * Watt * PetscPowRealInt(meter, -2) * PetscPowRealInt(Kelvin, -4); /*… in RDCreate()
1090 rd->MaterialEnergy = RDMaterialEnergy_Saha; in RDCreate()
1094 …(DMDACreate1d(comm, DM_BOUNDARY_NONE, 20, sizeof(RDNode) / sizeof(PetscScalar), 1, NULL, &rd->da)); in RDCreate()
1095 PetscCall(DMSetFromOptions(rd->da)); in RDCreate()
1096 PetscCall(DMSetUp(rd->da)); in RDCreate()
1097 PetscCall(DMDASetFieldName(rd->da, 0, "E")); in RDCreate()
1098 PetscCall(DMDASetFieldName(rd->da, 1, "T")); in RDCreate()
1099 PetscCall(DMDASetUniformCoordinates(rd->da, 0., 1., 0., 0., 0., 0.)); in RDCreate()
1101 *inrd = rd; in RDCreate()
1107 RD rd; in main() local
1117 PetscCall(RDCreate(PETSC_COMM_WORLD, &rd)); in main()
1118 PetscCall(DMCreateGlobalVector(rd->da, &X)); in main()
1119 PetscCall(DMSetMatType(rd->da, MATAIJ)); in main()
1120 PetscCall(DMCreateMatrix(rd->da, &B)); in main()
1121 PetscCall(RDInitialState(rd, X)); in main()
1126 PetscCall(TSSetDM(ts, rd->da)); in main()
1127 switch (rd->discretization) { in main()
1129 PetscCall(TSSetIFunction(ts, NULL, RDIFunction_FD, rd)); in main()
1130 if (rd->jacobian == JACOBIAN_ANALYTIC) PetscCall(TSSetIJacobian(ts, B, B, RDIJacobian_FD, rd)); in main()
1133 PetscCall(TSSetIFunction(ts, NULL, RDIFunction_FE, rd)); in main()
1134 if (rd->jacobian == JACOBIAN_ANALYTIC) PetscCall(TSSetIJacobian(ts, B, B, RDIJacobian_FE, rd)); in main()
1137 PetscCall(TSSetMaxTime(ts, rd->final_time)); in main()
1144 switch (rd->jacobian) { in main()
1157 if (rd->test_diff) PetscCall(RDTestDifferentiation(rd)); in main()
1162 if (rd->view_draw) PetscCall(RDView(rd, X, PETSC_VIEWER_DRAW_WORLD)); in main()
1163 if (rd->view_binary[0]) { in main()
1165 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, rd->view_binary, FILE_MODE_WRITE, &viewer)); in main()
1166 PetscCall(RDView(rd, X, viewer)); in main()
1171 PetscCall(RDDestroy(&rd)); in main()