Actual source code: ex2.c
1: static char help[] = "Runaway electron model with Landau collision operator\n\n";
3: #include <petscdmplex.h>
4: #include <petsclandau.h>
5: #include <petscts.h>
6: #include <petscds.h>
7: #include <petscdmcomposite.h>
9: /* data for runaway electron model */
10: typedef struct REctx_struct {
11: PetscErrorCode (*test)(TS, Vec, PetscInt, PetscReal, PetscBool, LandauCtx *, struct REctx_struct *);
12: PetscErrorCode (*impuritySrcRate)(PetscReal, PetscReal *, LandauCtx*);
13: PetscErrorCode (*E)(Vec, Vec, PetscInt, PetscReal, LandauCtx*, PetscReal *);
14: PetscReal T_cold; /* temperature of newly ionized electrons and impurity ions */
15: PetscReal ion_potential; /* ionization potential of impurity */
16: PetscReal Ne_ion; /* effective number of electrons shed in ioization of impurity */
17: PetscReal Ez_initial;
18: PetscReal L; /* inductance */
19: Vec X_0;
20: PetscInt imp_idx; /* index for impurity ionizing sink */
21: PetscReal pulse_start;
22: PetscReal pulse_width;
23: PetscReal pulse_rate;
24: PetscReal current_rate;
25: PetscInt plotIdx;
26: PetscInt plotStep;
27: PetscInt idx; /* cache */
28: PetscReal j; /* cache */
29: PetscReal plotDt;
30: PetscBool plotting;
31: PetscBool use_spitzer_eta;
32: PetscInt print_period;
33: } REctx;
35: static const PetscReal kev_joul = 6.241506479963235e+15; /* 1/1000e */
37: #define RE_CUT 3.
38: /* < v, u_re * v * q > */
39: static void f0_j_re(PetscInt dim, PetscInt Nf, PetscInt NfAux,
40: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
41: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
42: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
43: {
44: PetscReal n_e = PetscRealPart(u[0]);
45: if (dim==2) {
46: if (x[1] > RE_CUT || x[1] < -RE_CUT) { /* simply a cutoff for REs. v_|| > 3 v(T_e) */
47: *f0 = n_e * 2.*PETSC_PI*x[0] * x[1] * constants[0]; /* n * r * v_|| * q */
48: } else {
49: *f0 = 0;
50: }
51: } else {
52: if (x[2] > RE_CUT || x[2] < -RE_CUT) { /* simply a cutoff for REs. v_|| > 3 v(T_e) */
53: *f0 = n_e * x[2] * constants[0];
54: } else {
55: *f0 = 0;
56: }
57: }
58: }
60: /* sum < v, u*v*q > */
61: static void f0_jz_sum(PetscInt dim, PetscInt Nf, PetscInt NfAux,
62: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
63: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
64: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar q[], PetscScalar *f0)
65: {
66: PetscInt ii;
67: f0[0] = 0;
68: if (dim==2) {
69: for (ii=0;ii<Nf;ii++) f0[0] += u[ii] * 2.*PETSC_PI*x[0] * x[1] * q[ii]; /* n * r * v_|| * q * v_0 */
70: } else {
71: for (ii=0;ii<Nf;ii++) f0[0] += u[ii] * x[2] * q[ii]; /* n * v_|| * q * v_0 */
72: }
73: }
75: /* < v, n_e > */
76: static void f0_n(PetscInt dim, PetscInt Nf, PetscInt NfAux,
77: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
78: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
79: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
80: {
81: PetscInt ii = (PetscInt)PetscRealPart(constants[0]);
82: if (dim==2) f0[0] = 2.*PETSC_PI*x[0]*u[ii];
83: else {
84: f0[0] = u[ii];
85: }
86: }
88: /* < v, n_e v_|| > */
89: static void f0_vz(PetscInt dim, PetscInt Nf, PetscInt NfAux,
90: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
91: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
92: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
93: {
94: PetscInt ii = (PetscInt)PetscRealPart(constants[0]);
95: if (dim==2) f0[0] = u[ii] * 2.*PETSC_PI*x[0] * x[1]; /* n r v_|| */
96: else {
97: f0[0] = u[ii] * x[2]; /* n v_|| */
98: }
99: }
101: /* < v, n_e (v-shift) > */
102: static void f0_ve_shift(PetscInt dim, PetscInt Nf, PetscInt NfAux,
103: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
104: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
105: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
106: {
107: PetscReal vz = numConstants>0 ? PetscRealPart(constants[0]) : 0;
108: if (dim==2) *f0 = u[0] * 2.*PETSC_PI*x[0] * PetscSqrtReal(x[0]*x[0] + (x[1]-vz)*(x[1]-vz)); /* n r v */
109: else {
110: *f0 = u[0] * PetscSqrtReal(x[0]*x[0] + x[1]*x[1] + (x[2]-vz)*(x[2]-vz)); /* n v */
111: }
112: }
114: /* CalculateE - Calculate the electric field */
115: /* T -- Electron temperature */
116: /* n -- Electron density */
117: /* lnLambda -- */
118: /* eps0 -- */
119: /* E -- output E, input \hat E */
120: static PetscReal CalculateE(PetscReal Tev, PetscReal n, PetscReal lnLambda, PetscReal eps0, PetscReal *E)
121: {
122: PetscReal c,e,m;
125: c = 299792458.0;
126: e = 1.602176e-19;
127: m = 9.10938e-31;
128: if (1) {
129: double Ec, Ehat = *E, betath = PetscSqrtReal(2*Tev*e/(m*c*c)), j0 = Ehat * 7/(PetscSqrtReal(2)*2) * PetscPowReal(betath,3) * n * e * c;
130: Ec = n*lnLambda*PetscPowReal(e,3) / (4*PETSC_PI*PetscPowReal(eps0,2)*m*c*c);
131: *E = Ec;
132: PetscPrintf(PETSC_COMM_WORLD, "CalculateE j0=%g Ec = %g\n",j0,Ec);
133: } else {
134: PetscReal Ed, vth;
135: vth = PetscSqrtReal(8*Tev*e/(m*PETSC_PI));
136: Ed = n*lnLambda*PetscPowReal(e,3) / (4*PETSC_PI*PetscPowReal(eps0,2)*m*vth*vth);
137: *E = Ed;
138: }
139: return(0);
140: }
142: static PetscReal Spitzer(PetscReal m_e, PetscReal e, PetscReal Z, PetscReal epsilon0, PetscReal lnLam, PetscReal kTe_joules)
143: {
144: PetscReal Fz = (1+1.198*Z+0.222*Z*Z)/(1+2.966*Z+0.753*Z*Z), eta;
145: eta = Fz*4./3.*PetscSqrtReal(2.*PETSC_PI)*Z*PetscSqrtReal(m_e)*PetscSqr(e)*lnLam*PetscPowReal(4*PETSC_PI*epsilon0,-2.)*PetscPowReal(kTe_joules,-1.5);
146: return eta;
147: }
149: /* */
150: static PetscErrorCode testNone(TS ts, Vec X, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
151: {
153: return(0);
154: }
156: /* */
157: static PetscErrorCode testSpitzer(TS ts, Vec X, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
158: {
159: PetscErrorCode ierr;
160: PetscInt ii;
161: PetscDS prob;
162: static PetscReal old_ratio = 1e10;
163: TSConvergedReason reason;
164: PetscReal J,J_re,spit_eta,Te_kev=0,E,ratio,Z,n_e,v,v2;
165: PetscScalar user[2] = {0.,ctx->charges[0]}, q[LANDAU_MAX_SPECIES],tt[LANDAU_MAX_SPECIES],vz;
166: PetscReal dt;
167: DM pack, plexe = ctx->plex[0], plexi = (ctx->num_grids==1) ? NULL : ctx->plex[1];
168: Vec XsubArray[LANDAU_MAX_GRIDS];
171: if (ctx->num_species!=2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "ctx->num_species %D != 2",ctx->num_species);
172: VecGetDM(X, &pack);
173: if (!pack) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "no DM");
174: DMCompositeGetAccessArray(pack, X, ctx->num_grids, NULL, XsubArray); // read only
175: TSGetTimeStep(ts,&dt);
176: /* get current for each grid */
177: for (ii=0;ii<ctx->num_species;ii++) q[ii] = ctx->charges[ii];
178: DMGetDS(plexe, &prob);
179: PetscDSSetConstants(prob, 2, &q[0]);
180: PetscDSSetObjective(prob, 0, &f0_jz_sum);
181: DMPlexComputeIntegralFEM(plexe,XsubArray[0],tt,NULL);
182: J = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]);
183: if (plexi) { // add ions
184: DMGetDS(plexi, &prob);
185: PetscDSSetConstants(prob, 1, &q[1]);
186: PetscDSSetObjective(prob, 0, &f0_jz_sum);
187: DMPlexComputeIntegralFEM(plexi,XsubArray[1],tt,NULL);
188: J += -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]);
189: }
190: PetscPrintf(ctx->comm, "testSpitzer J = %10.3e\n",J);
191: /* get N_e */
192: DMGetDS(plexe, &prob);
193: PetscDSSetConstants(prob, 1, user);
194: PetscDSSetObjective(prob, 0, &f0_n);
195: DMPlexComputeIntegralFEM(plexe,XsubArray[0],tt,NULL);
196: n_e = PetscRealPart(tt[0])*ctx->n_0;
197: /* Z */
198: Z = -ctx->charges[1]/ctx->charges[0];
199: /* remove drift */
200: if (0) {
201: user[0] = 0; // electrons
202: DMGetDS(plexe, &prob);
203: PetscDSSetConstants(prob, 1, user);
204: PetscDSSetObjective(prob, 0, &f0_vz);
205: DMPlexComputeIntegralFEM(plexe,XsubArray[0],tt,NULL);
206: vz = ctx->n_0*PetscRealPart(tt[0])/n_e; /* non-dimensional */
207: } else vz = 0;
208: /* thermal velocity */
209: DMGetDS(plexe, &prob);
210: PetscDSSetConstants(prob, 1, &vz);
211: PetscDSSetObjective(prob, 0, &f0_ve_shift);
212: DMPlexComputeIntegralFEM(plexe,XsubArray[0],tt,NULL);
213: v = ctx->n_0*ctx->v_0*PetscRealPart(tt[0])/n_e; /* remove number density to get velocity */
214: v2 = PetscSqr(v); /* use real space: m^2 / s^2 */
215: Te_kev = (v2*ctx->masses[0]*PETSC_PI/8)*kev_joul; /* temperature in kev */
216: spit_eta = Spitzer(ctx->masses[0],-ctx->charges[0],Z,ctx->epsilon0,ctx->lnLam,Te_kev/kev_joul); /* kev --> J (kT) */
217: if (0) {
218: DMGetDS(plexe, &prob);
219: PetscDSSetConstants(prob, 1, q);
220: PetscDSSetObjective(prob, 0, &f0_j_re);
221: DMPlexComputeIntegralFEM(plexe,XsubArray[0],tt,NULL);
222: } else tt[0] = 0;
223: J_re = -ctx->n_0*ctx->v_0*PetscRealPart(tt[0]);
224: DMCompositeRestoreAccessArray(pack, X, ctx->num_grids, NULL, XsubArray); // read only
226: if (rectx->use_spitzer_eta) {
227: E = ctx->Ez = spit_eta*(rectx->j-J_re);
228: } else {
229: E = ctx->Ez; /* keep real E */
230: rectx->j = J; /* cache */
231: }
233: ratio = E/J/spit_eta;
234: if (stepi>10 && !rectx->use_spitzer_eta && (
235: //(old_ratio-ratio < 1.e-3 && ratio > 0.99 && ratio < 1.01) ||
236: //(old_ratio-ratio < 1.e-4 && ratio > 0.98 && ratio < 1.02) ||
237: (old_ratio-ratio < 1.e-6))) {
238: rectx->pulse_start = time + 0.98*dt;
239: rectx->use_spitzer_eta = PETSC_TRUE;
240: }
241: TSGetConvergedReason(ts,&reason);
242: TSGetConvergedReason(ts,&reason);
243: if ((rectx->plotting) || stepi == 0 || reason || rectx->pulse_start == time + 0.98*dt) {
244: PetscPrintf(ctx->comm, "testSpitzer: %4D) time=%11.4e n_e= %10.3e E= %10.3e J= %10.3e J_re= %10.3e %.3g%% Te_kev= %10.3e Z_eff=%g E/J to eta ratio= %g (diff=%g) %s %s\n",stepi,time,n_e/ctx->n_0,ctx->Ez,J,J_re,100*J_re/J, Te_kev,Z,ratio,old_ratio-ratio, rectx->use_spitzer_eta ? "using Spitzer eta*J E" : "constant E",rectx->pulse_start != time + 0.98*dt ? "normal" : "transition");
245: if (rectx->pulse_start == time + 0.98*dt) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Spitzer complet ratio=%g",ratio);
246: }
247: old_ratio = ratio;
248: return(0);
249: }
251: static const double ppp = 2;
252: static void f0_0_diff_lp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
253: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
254: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
255: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
256: {
257: LandauCtx *ctx = (LandauCtx *)constants;
258: REctx *rectx = (REctx*)ctx->data;
259: PetscInt ii = rectx->idx, i;
260: const PetscReal kT_m = ctx->k*ctx->thermal_temps[ii]/ctx->masses[ii]; /* kT/m */
261: const PetscReal n = ctx->n[ii];
262: PetscReal diff, f_maxwell, v2 = 0, theta = 2*kT_m/(ctx->v_0*ctx->v_0); /* theta = 2kT/mc^2 */
263: for (i = 0; i < dim; ++i) v2 += x[i]*x[i];
264: f_maxwell = n*PetscPowReal(PETSC_PI*theta,-1.5)*(PetscExpReal(-v2/theta));
265: diff = 2.*PETSC_PI*x[0]*(PetscRealPart(u[ii]) - f_maxwell);
266: f0[0] = PetscPowReal(diff,ppp);
267: }
268: static void f0_0_maxwellian_lp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
269: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
270: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
271: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar *f0)
272: {
273: LandauCtx *ctx = (LandauCtx *)constants;
274: REctx *rectx = (REctx*)ctx->data;
275: PetscInt ii = rectx->idx, i;
276: const PetscReal kT_m = ctx->k*ctx->thermal_temps[ii]/ctx->masses[ii]; /* kT/m */
277: const PetscReal n = ctx->n[ii];
278: PetscReal f_maxwell, v2 = 0, theta = 2*kT_m/(ctx->v_0*ctx->v_0); /* theta = 2kT/mc^2 */
279: for (i = 0; i < dim; ++i) v2 += x[i]*x[i];
280: f_maxwell = 2.*PETSC_PI*x[0] * n*PetscPowReal(PETSC_PI*theta,-1.5)*(PetscExpReal(-v2/theta));
281: f0[0] = PetscPowReal(f_maxwell,ppp);
282: }
284: /* */
285: static PetscErrorCode testStable(TS ts, Vec X, PetscInt stepi, PetscReal time, PetscBool islast, LandauCtx *ctx, REctx *rectx)
286: {
287: PetscErrorCode ierr;
288: PetscDS prob;
289: Vec X2;
290: PetscReal ediff,idiff=0,lpm0,lpm1=1;
291: PetscScalar tt[LANDAU_MAX_SPECIES];
292: DM dm, plex = ctx->plex[0];
295: VecGetDM(X, &dm);
296: DMGetDS(plex, &prob);
297: VecDuplicate(X,&X2);
298: VecCopy(X,X2);
299: if (!rectx->X_0) {
300: VecDuplicate(X,&rectx->X_0);
301: VecCopy(X,rectx->X_0);
302: }
303: VecAXPY(X,-1.0,rectx->X_0);
304: PetscDSSetConstants(prob, sizeof(LandauCtx)/sizeof(PetscScalar), (PetscScalar*)ctx);
305: rectx->idx = 0;
306: PetscDSSetObjective(prob, 0, &f0_0_diff_lp);
307: DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
308: ediff = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
309: PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp);
310: DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
311: lpm0 = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
312: if (ctx->num_species>1) {
313: rectx->idx = 1;
314: PetscDSSetObjective(prob, 0, &f0_0_diff_lp);
315: DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
316: idiff = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
317: PetscDSSetObjective(prob, 0, &f0_0_maxwellian_lp);
318: DMPlexComputeIntegralFEM(plex,X2,tt,NULL);
319: lpm1 = PetscPowReal(PetscRealPart(tt[0]),1./ppp);
320: }
321: PetscPrintf(PETSC_COMM_WORLD, "%s %D) time=%10.3e n-%d norm electrons/max=%20.13e ions/max=%20.13e\n", "----",stepi,time,(int)ppp,ediff/lpm0,idiff/lpm1);
322: /* view */
323: VecCopy(X2,X);
324: VecDestroy(&X2);
325: if (islast) {
326: VecDestroy(&rectx->X_0);
327: rectx->X_0 = NULL;
328: }
329: return(0);
330: }
332: static PetscErrorCode EInduction(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
333: {
334: REctx *rectx = (REctx*)ctx->data;
335: PetscErrorCode ierr;
336: PetscInt ii;
337: DM dm,plex;
338: PetscScalar tt[LANDAU_MAX_SPECIES], qv0[LANDAU_MAX_SPECIES];
339: PetscReal dJ_dt;
340: PetscDS prob;
343: for (ii=0;ii<ctx->num_species;ii++) qv0[ii] = ctx->charges[ii]*ctx->v_0;
344: VecGetDM(X, &dm);
345: DMGetDS(dm, &prob);
346: DMConvert(dm, DMPLEX, &plex);
347: /* get d current / dt */
348: PetscDSSetConstants(prob, ctx->num_species, qv0);
349: PetscDSSetObjective(prob, 0, &f0_jz_sum);
350: if (!X_t) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "X_t");
351: DMPlexComputeIntegralFEM(plex,X_t,tt,NULL);
352: dJ_dt = -ctx->n_0*PetscRealPart(tt[0])/ctx->t_0;
353: /* E induction */
354: *a_E = -rectx->L*dJ_dt + rectx->Ez_initial;
355: DMDestroy(&plex);
356: return(0);
357: }
359: static PetscErrorCode EConstant(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
360: {
362: *a_E = ctx->Ez;
363: return(0);
364: }
366: static PetscErrorCode ENone(Vec X, Vec X_t, PetscInt step, PetscReal time, LandauCtx *ctx, PetscReal *a_E)
367: {
369: *a_E = 0;
370: return(0);
371: }
373: /* ------------------------------------------------------------------- */
374: /*
375: FormSource - Evaluates source terms F(t).
377: Input Parameters:
378: . ts - the TS context
379: . time -
380: . X_dummmy - input vector
381: . dummy - optional user-defined context, as set by SNESSetFunction()
383: Output Parameter:
384: . F - function vector
385: */
386: static PetscErrorCode FormSource(TS ts, PetscReal ftime, Vec X_dummmy, Vec F, void *dummy)
387: {
388: PetscReal new_imp_rate;
389: LandauCtx *ctx;
390: DM pack;
392: REctx *rectx;
395: TSGetDM(ts,&pack);
396: DMGetApplicationContext(pack, &ctx);
397: rectx = (REctx*)ctx->data;
398: /* check for impurities */
399: rectx->impuritySrcRate(ftime,&new_imp_rate,ctx);
400: if (new_imp_rate != 0) {
401: if (new_imp_rate != rectx->current_rate) {
402: PetscInt ii;
403: PetscReal dne_dt,dni_dt,tilda_ns[LANDAU_MAX_SPECIES],temps[LANDAU_MAX_SPECIES];
404: Vec globFarray[LANDAU_MAX_GRIDS];
405: rectx->current_rate = new_imp_rate;
406: for (ii=1;ii<LANDAU_MAX_SPECIES;ii++) tilda_ns[ii] = 0;
407: for (ii=1;ii<LANDAU_MAX_SPECIES;ii++) temps[ii] = 1;
408: dni_dt = new_imp_rate /* *ctx->t_0 */; /* fully ionized immediately, no normalize, stay in non-dim */
409: dne_dt = new_imp_rate*rectx->Ne_ion /* *ctx->t_0 */;
410: tilda_ns[0] = dne_dt; tilda_ns[rectx->imp_idx] = dni_dt;
411: temps[0] = rectx->T_cold; temps[rectx->imp_idx] = rectx->T_cold;
412: PetscInfo4(ctx->plex[0], "\tHave new_imp_rate= %10.3e time= %10.3e de/dt= %10.3e di/dt= %10.3e ***\n",new_imp_rate,ftime,dne_dt,dni_dt);
413: DMCompositeGetAccessArray(pack, F, ctx->num_grids, NULL, globFarray);
414: for (PetscInt grid=0 ; grid<ctx->num_grids ; grid++) {
415: /* add it */
416: LandauAddMaxwellians(ctx->plex[grid],globFarray[grid],ftime,temps,tilda_ns,grid,ctx);
417: VecViewFromOptions(globFarray[grid],NULL,"-vec_view_sources");
418: }
419: DMCompositeRestoreAccessArray(pack, F, ctx->num_grids, NULL, globFarray);
420: }
421: } else {
422: VecZeroEntries(F);
423: rectx->current_rate = 0;
424: }
425: return(0);
426: }
427: PetscErrorCode Monitor(TS ts, PetscInt stepi, PetscReal time, Vec X, void *actx)
428: {
429: LandauCtx *ctx = (LandauCtx*) actx; /* user-defined application context */
430: REctx *rectx = (REctx*)ctx->data;
431: DM pack;
432: Vec globXArray[LANDAU_MAX_GRIDS];
433: TSConvergedReason reason;
434: PetscErrorCode ierr;
436: VecGetDM(X, &pack);
437: DMCompositeGetAccessArray(pack, X, ctx->num_grids, NULL, globXArray);
438: if (stepi > rectx->plotStep && rectx->plotting) {
439: rectx->plotting = PETSC_FALSE; /* was doing diagnostics, now done */
440: rectx->plotIdx++;
441: }
442: /* view */
443: TSGetConvergedReason(ts,&reason);
444: if (time/rectx->plotDt >= (PetscReal)rectx->plotIdx || reason) {
445: if ((reason || stepi==0 || rectx->plotIdx%rectx->print_period==0) && ctx->verbose > 0) {
446: /* print norms */
447: LandauPrintNorms(X, stepi);
448: }
449: if (!rectx->plotting) { /* first step of possible backtracks */
450: rectx->plotting = PETSC_TRUE;
451: /* diagnostics + change E field with Sptizer (not just a monitor) */
452: rectx->test(ts,X,stepi,time,reason ? PETSC_TRUE : PETSC_FALSE, ctx, rectx);
453: } else {
454: PetscPrintf(PETSC_COMM_WORLD, "\t\t ERROR SKIP test spit ------\n");
455: rectx->plotting = PETSC_TRUE;
456: }
457: PetscObjectSetName((PetscObject) globXArray[0], ctx->plex[1] ? "ue" : "ux");
458: if (ctx->plex[1]) {
459: PetscObjectSetName((PetscObject) globXArray[1], "ui");
460: }
461: /* view, overwrite step when back tracked */
462: DMSetOutputSequenceNumber(pack, rectx->plotIdx, time*ctx->t_0);
463: VecViewFromOptions(globXArray[0],NULL,"-vec_view_e");
464: if (ctx->plex[1]) {
465: VecViewFromOptions(globXArray[1],NULL,"-vec_view_i");
466: }
467: rectx->plotStep = stepi;
468: } else {
469: if (rectx->plotting) PetscPrintf(PETSC_COMM_WORLD," ERROR rectx->plotting=%D step %D\n",rectx->plotting,stepi);
470: /* diagnostics + change E field with Sptizer (not just a monitor) - can we lag this? */
471: rectx->test(ts,X,stepi,time,reason ? PETSC_TRUE : PETSC_FALSE, ctx, rectx);
472: }
473: DMCompositeRestoreAccessArray(pack, X, ctx->num_grids, NULL, globXArray);
474: /* parallel check */
475: if (reason && ctx->verbose > 0) {
476: PetscReal val,rval;
477: PetscMPIInt rank;
478: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
479: TSGetSolution(ts, &X);
480: VecNorm(X,NORM_2,&val);
481: MPIU_Allreduce(&val,&rval,1,MPIU_REAL,MPIU_MAX,PETSC_COMM_WORLD);
482: if (rval != val) {
483: PetscPrintf(PETSC_COMM_SELF, " ***** [%D] ERROR max |x| = %22.15e, my |x| = %22.15e diff=%e\n",rank,rval,val,rval-val);
484: } else {
485: PetscPrintf(PETSC_COMM_WORLD, "[%D] parallel consistency check OK\n",rank);
486: }
487: }
488: rectx->idx = 0;
489: return(0);
490: }
492: PetscErrorCode PreStep(TS ts)
493: {
495: LandauCtx *ctx;
496: REctx *rectx;
497: DM dm;
498: PetscInt stepi;
499: PetscReal time;
500: Vec X;
503: /* not used */
504: TSGetDM(ts,&dm);
505: TSGetTime(ts,&time);
506: TSGetSolution(ts,&X);
507: DMGetApplicationContext(dm, &ctx);
508: rectx = (REctx*)ctx->data;
509: TSGetStepNumber(ts, &stepi);
510: /* update E */
511: rectx->E(X, NULL, stepi, time, ctx, &ctx->Ez);
512: return(0);
513: }
515: /* model for source of non-ionized impurities, profile provided by model, in du/dt form in normalized units (tricky because n_0 is normalized with electrons) */
516: static PetscErrorCode stepSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
517: {
518: REctx *rectx = (REctx*)ctx->data;
521: if (time >= rectx->pulse_start) *rho = rectx->pulse_rate;
522: else *rho = 0.;
523: return(0);
524: }
525: static PetscErrorCode zeroSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
526: {
528: *rho = 0.;
529: return(0);
530: }
531: static PetscErrorCode pulseSrc(PetscReal time, PetscReal *rho, LandauCtx *ctx)
532: {
533: REctx *rectx = (REctx*)ctx->data;
536: if (rectx->pulse_start == PETSC_MAX_REAL) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"'-ex2_pulse_start_time X' must be used with '-ex2_impurity_source_type pulse'");
537: if (time < rectx->pulse_start || time > rectx->pulse_start + 3*rectx->pulse_width) *rho = 0;
538: /* else if (0) { */
539: /* double t = time - rectx->pulse_start, start = rectx->pulse_width, stop = 2*rectx->pulse_width, cycle = 3*rectx->pulse_width, steep = 5, xi = 0.75 - (stop - start)/(2* cycle); */
540: /* *rho = rectx->pulse_rate * (cycle / (stop - start)) / (1 + PetscExpReal(steep*(PetscSinReal(2*PETSC_PI*((t - start)/cycle + xi)) - PetscSinReal(2*PETSC_PI*xi)))); */
541: /* } else if (0) { */
542: /* double x = 2*(time - rectx->pulse_start)/(3*rectx->pulse_width) - 1; */
543: /* if (x==1 || x==-1) *rho = 0; */
544: /* else *rho = rectx->pulse_rate * PetscExpReal(-1/(1-x*x)); */
545: /* } */
546: else {
547: double x = PetscSinReal((time-rectx->pulse_start)/(3*rectx->pulse_width)*2*PETSC_PI - PETSC_PI/2) + 1; /* 0:2, integrates to 1.0 */
548: *rho = rectx->pulse_rate * x / (3*rectx->pulse_width);
549: if (!rectx->use_spitzer_eta) rectx->use_spitzer_eta = PETSC_TRUE; /* use it next time */
550: }
551: return(0);
552: }
556: static PetscErrorCode ProcessREOptions(REctx *rectx, const LandauCtx *ctx, DM dm, const char prefix[])
557: {
558: PetscErrorCode ierr;
559: PetscFunctionList plist = NULL, testlist = NULL, elist = NULL;
560: char pname[256],testname[256],ename[256];
561: DM dm_dummy;
562: PetscBool Connor_E = PETSC_FALSE;
565: DMCreate(PETSC_COMM_WORLD,&dm_dummy);
566: rectx->Ne_ion = 1; /* number of electrons given up by impurity ion */
567: rectx->T_cold = .005; /* kev */
568: rectx->ion_potential = 15; /* ev */
569: rectx->L = 2;
570: rectx->X_0 = NULL;
571: rectx->imp_idx = ctx->num_species - 1; /* default ionized impurity as last one */
572: rectx->pulse_start = PETSC_MAX_REAL;
573: rectx->pulse_width = 1;
574: rectx->plotStep = PETSC_MAX_INT;
575: rectx->pulse_rate = 1.e-1;
576: rectx->current_rate = 0;
577: rectx->plotIdx = 0;
578: rectx->j = 0;
579: rectx->plotDt = 1.0;
580: rectx->plotting = PETSC_FALSE;
581: rectx->use_spitzer_eta = PETSC_FALSE;
582: rectx->idx = 0;
583: rectx->print_period = 10;
584: /* Register the available impurity sources */
585: PetscFunctionListAdd(&plist,"step",&stepSrc);
586: PetscFunctionListAdd(&plist,"none",&zeroSrc);
587: PetscFunctionListAdd(&plist,"pulse",&pulseSrc);
588: PetscStrcpy(pname,"none");
589: PetscFunctionListAdd(&testlist,"none",&testNone);
590: PetscFunctionListAdd(&testlist,"spitzer",&testSpitzer);
591: PetscFunctionListAdd(&testlist,"stable",&testStable);
592: PetscStrcpy(testname,"none");
593: PetscFunctionListAdd(&elist,"none",&ENone);
594: PetscFunctionListAdd(&elist,"induction",&EInduction);
595: PetscFunctionListAdd(&elist,"constant",&EConstant);
596: PetscStrcpy(ename,"constant");
598: PetscOptionsBegin(PETSC_COMM_SELF, prefix, "Options for Runaway/seed electron model", "none");
599: PetscOptionsReal("-ex2_plot_dt", "Plotting interval", "xgc_dmplex.c", rectx->plotDt, &rectx->plotDt, NULL);
600: if (rectx->plotDt < 0) rectx->plotDt = 1e30;
601: if (rectx->plotDt == 0) rectx->plotDt = 1e-30;
602: PetscOptionsInt("-ex2_print_period", "Plotting interval", "xgc_dmplex.c", rectx->print_period, &rectx->print_period, NULL);
603: PetscOptionsFList("-ex2_impurity_source_type","Name of impurity source to run","",plist,pname,pname,sizeof(pname),NULL);
604: PetscOptionsFList("-ex2_test_type","Name of test to run","",testlist,testname,testname,sizeof(testname),NULL);
605: PetscOptionsInt("-ex2_impurity_index", "index of sink for impurities", "none", rectx->imp_idx, &rectx->imp_idx, NULL);
606: if ((rectx->imp_idx >= ctx->num_species || rectx->imp_idx < 1) && ctx->num_species > 1) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"index of sink for impurities ions is out of range (%D), must be > 0 && < NS",rectx->imp_idx);
607: PetscOptionsFList("-ex2_e_field_type","Electric field type","",elist,ename,ename,sizeof(ename),NULL);
608: rectx->Ne_ion = -ctx->charges[rectx->imp_idx]/ctx->charges[0];
609: PetscOptionsReal("-ex2_t_cold","Temperature of cold electron and ions after ionization in keV","none",rectx->T_cold,&rectx->T_cold, NULL);
610: PetscOptionsReal("-ex2_pulse_start_time","Time at which pulse happens for 'pulse' source","none",rectx->pulse_start,&rectx->pulse_start, NULL);
611: PetscOptionsReal("-ex2_pulse_width_time","Width of pulse 'pulse' source","none",rectx->pulse_width,&rectx->pulse_width, NULL);
612: PetscOptionsReal("-ex2_pulse_rate","Number density of pulse for 'pulse' source","none",rectx->pulse_rate,&rectx->pulse_rate, NULL);
613: rectx->T_cold *= 1.16e7; /* convert to Kelvin */
614: PetscOptionsReal("-ex2_ion_potential","Potential to ionize impurity (should be array) in ev","none",rectx->ion_potential,&rectx->ion_potential, NULL);
615: PetscOptionsReal("-ex2_inductance","Inductance E feild","none",rectx->L,&rectx->L, NULL);
616: PetscOptionsBool("-ex2_connor_e_field_units","Scale Ex but Connor-Hastie E_c","none",Connor_E,&Connor_E, NULL);
617: PetscInfo5(dm_dummy, "Num electrons from ions=%g, T_cold=%10.3e, ion potential=%10.3e, E_z=%10.3e v_0=%10.3e\n",rectx->Ne_ion,rectx->T_cold,rectx->ion_potential,ctx->Ez,ctx->v_0);
618: PetscOptionsEnd();
619: /* get impurity source rate function */
620: PetscFunctionListFind(plist,pname,&rectx->impuritySrcRate);
621: if (!rectx->impuritySrcRate) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No impurity source function found '%s'",pname);
622: PetscFunctionListFind(testlist,testname,&rectx->test);
623: if (!rectx->test) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No test found '%s'",testname);
624: PetscFunctionListFind(elist,ename,&rectx->E);
625: if (!rectx->E) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"No E field function found '%s'",ename);
626: PetscFunctionListDestroy(&plist);
627: PetscFunctionListDestroy(&testlist);
628: PetscFunctionListDestroy(&elist);
630: /* convert E from Connor-Hastie E_c units to real if doing Spitzer E */
631: if (Connor_E) {
632: PetscReal E = ctx->Ez, Tev = ctx->thermal_temps[0]*8.621738e-5, n = ctx->n_0*ctx->n[0];
633: CalculateE(Tev, n, ctx->lnLam, ctx->epsilon0, &E);
634: ((LandauCtx*)ctx)->Ez *= E;
635: }
636: DMDestroy(&dm_dummy);
637: return(0);
638: }
640: int main(int argc, char **argv)
641: {
642: DM pack;
643: Vec X,XsubArray[LANDAU_MAX_GRIDS];
645: PetscInt dim = 2;
646: TS ts;
647: Mat J;
648: PetscDS prob;
649: LandauCtx *ctx;
650: REctx *rectx;
651: #if defined PETSC_USE_LOG
652: PetscLogStage stage;
653: #endif
654: PetscMPIInt rank;
655: PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
656: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
657: if (rank) { /* turn off output stuff for duplicate runs */
658: PetscOptionsClearValue(NULL,"-dm_view_e");
659: PetscOptionsClearValue(NULL,"-vec_view_e");
660: PetscOptionsClearValue(NULL,"-dm_view_diff_e");
661: PetscOptionsClearValue(NULL,"-vec_view_diff_e");
662: PetscOptionsClearValue(NULL,"-dm_view_sources_e");
663: PetscOptionsClearValue(NULL,"-vec_view_0_e");
664: PetscOptionsClearValue(NULL,"-dm_view_0_e");
665: PetscOptionsClearValue(NULL,"-vec_view_sources_e");
666: PetscOptionsClearValue(NULL,"-info"); /* this does not work */
667: PetscOptionsClearValue(NULL,"-dm_view_i");
668: PetscOptionsClearValue(NULL,"-vec_view_i");
669: PetscOptionsClearValue(NULL,"-dm_view_diff_i");
670: PetscOptionsClearValue(NULL,"-vec_view_diff_i");
671: PetscOptionsClearValue(NULL,"-dm_view_sources_i");
672: PetscOptionsClearValue(NULL,"-vec_view_sources_i");
673: PetscOptionsClearValue(NULL,"-dm_view_0_i");
674: PetscOptionsClearValue(NULL,"-vec_view_0_i");
675: }
676: PetscOptionsGetInt(NULL,NULL, "-dim", &dim, NULL);
677: /* Create a mesh */
678: LandauCreateVelocitySpace(PETSC_COMM_WORLD, dim, "", &X, &J, &pack);
679: PetscObjectSetName((PetscObject)J, "Jacobian");
680: PetscObjectSetName((PetscObject)X, "f");
681: LandauCreateMassMatrix(pack, NULL);
682: DMGetApplicationContext(pack, &ctx);
683: DMSetUp(pack);
684: DMGetDS(pack, &prob);
685: DMCompositeGetAccessArray(pack, X, ctx->num_grids, NULL, XsubArray); // read only
686: PetscObjectSetName((PetscObject) XsubArray[0], ctx->plex[1] ? "ue" : "ux");
687: if (ctx->plex[1]) {
688: PetscObjectSetName((PetscObject) XsubArray[1], "ui");
689: }
690: DMViewFromOptions(ctx->plex[0],NULL,"-dm_view_e");
691: DMViewFromOptions(ctx->plex[0],NULL,"-dm_view_0_e");
692: VecViewFromOptions(XsubArray[0],NULL,"-vec_view_0_e"); // inital condition (monitor plots after step)
693: if (ctx->plex[1]) {
694: DMViewFromOptions(ctx->plex[1],NULL,"-dm_view_i");
695: DMViewFromOptions(ctx->plex[1],NULL,"-dm_view_0_i");
696: VecViewFromOptions(XsubArray[1],NULL,"-vec_view_0_i"); // inital condition (monitor plots after step)
697: }
698: DMCompositeRestoreAccessArray(pack, X, ctx->num_grids, NULL, XsubArray); // read only
699: /* context */
700: PetscNew(&rectx);
701: ctx->data = rectx;
702: ProcessREOptions(rectx,ctx,pack,"");
703: DMSetOutputSequenceNumber(pack, 0, 0.0);
704: /* Create timestepping solver context */
705: TSCreate(PETSC_COMM_SELF,&ts);
706: TSSetDM(ts,pack);
707: TSSetIFunction(ts,NULL,LandauIFunction,NULL);
708: TSSetIJacobian(ts,J,J,LandauIJacobian,NULL);
709: TSSetRHSFunction(ts,NULL,FormSource,NULL);
710: TSSetFromOptions(ts);
711: TSSetSolution(ts,X);
712: TSSetApplicationContext(ts, ctx);
713: TSMonitorSet(ts,Monitor,ctx,NULL);
714: TSSetPreStep(ts,PreStep);
715: rectx->Ez_initial = ctx->Ez; /* cache for induction caclulation - applied E field */
716: if (1) { /* warm up an test just LandauIJacobian */
717: Vec vec;
718: PetscInt nsteps;
719: PetscReal dt;
720: PetscLogStageRegister("Warmup", &stage);
721: PetscLogStagePush(stage);
722: VecDuplicate(X,&vec);
723: VecCopy(X,vec);
724: TSGetMaxSteps(ts,&nsteps);
725: TSGetTimeStep(ts,&dt);
726: TSSetMaxSteps(ts,1);
727: TSSolve(ts,X);
728: TSSetMaxSteps(ts,nsteps);
729: TSSetStepNumber(ts,0);
730: TSSetTime(ts,0);
731: TSSetTimeStep(ts,dt);
732: rectx->plotIdx = 0;
733: rectx->plotting = PETSC_FALSE;
734: PetscLogStagePop();
735: VecCopy(vec,X);
736: VecDestroy(&vec);
737: ctx->aux_bool = PETSC_FALSE; // flag for not a clean Jacobian
738: }
739: /* go */
740: PetscLogStageRegister("Solve", &stage);
741: //TSSetSolution(ts,X);
742: MPI_Barrier(MPI_COMM_WORLD);
743: PetscLogStagePush(stage);
744: TSSolve(ts,X);
745: PetscLogStagePop();
746: /* clean up */
747: LandauDestroyVelocitySpace(&pack);
748: TSDestroy(&ts);
749: VecDestroy(&X);
750: PetscFree(rectx);
751: PetscFinalize();
752: return ierr;
753: }
755: /*TEST
756: test:
757: suffix: 0
758: requires: p4est !complex double
759: args: -dm_landau_num_species_grid 1,1 -dm_landau_Ez 0 -petscspace_degree 3 -petscspace_poly_tensor 1 -dm_landau_type p4est -info :dm,tsadapt -dm_landau_ion_masses 2 -dm_landau_ion_charges 1 -dm_landau_thermal_temps 5,5 -dm_landau_n 2,2 -dm_landau_n_0 5e19 -ts_monitor -snes_rtol 1.e-10 -snes_stol 1.e-14 -snes_monitor -snes_converged_reason -snes_max_it 10 -ts_type arkimex -ts_arkimex_type 1bee -ts_max_snes_failures -1 -ts_rtol 1e-3 -ts_dt 1.e-1 -ts_max_time 1 -ts_adapt_clip .5,1.25 -ts_max_steps 2 -ts_adapt_scale_solve_failed 0.75 -ts_adapt_time_step_increase_delay 5 -pc_type lu -ksp_type preonly -dm_landau_amr_levels_max 2,2 -ex2_impurity_source_type pulse -ex2_pulse_start_time 1e-1 -ex2_pulse_width_time 10 -ex2_pulse_rate 1e-2 -ex2_t_cold .05 -ex2_plot_dt 1e-1 -dm_refine 1 -dm_landau_gpu_assembly true -dm_landau_device_type cpu
761: test:
762: suffix: kokkos
763: requires: p4est !complex double kokkos_kernels
764: args: -dm_landau_num_species_grid 1,1 -dm_landau_Ez 0 -petscspace_degree 3 -petscspace_poly_tensor 1 -dm_landau_type p4est -info :tsadapt -dm_landau_ion_masses 2 -dm_landau_ion_charges 1 -dm_landau_thermal_temps 5,5 -dm_landau_n 2,2 -dm_landau_n_0 5e19 -ts_monitor -snes_rtol 1.e-10 -snes_stol 1.e-14 -snes_monitor -snes_converged_reason -snes_max_it 10 -ts_type arkimex -ts_arkimex_type 1bee -ts_max_snes_failures -1 -ts_rtol 1e-3 -ts_dt 1.e-1 -ts_max_time 1 -ts_adapt_clip .5,1.25 -ts_max_steps 2 -ts_adapt_scale_solve_failed 0.75 -ts_adapt_time_step_increase_delay 5 -pc_type lu -ksp_type preonly -dm_landau_amr_levels_max 3,2 -ex2_impurity_source_type pulse -ex2_pulse_start_time 1e-1 -ex2_pulse_width_time 10 -ex2_pulse_rate 1e-2 -ex2_t_cold .05 -ex2_plot_dt 1e-1 -dm_refine 1 -dm_landau_device_type kokkos -dm_landau_gpu_assembly true -dm_mat_type aijkokkos -dm_vec_type kokkos
766: test:
767: suffix: cuda
768: requires: p4est !complex double cuda
769: args: -dm_landau_num_species_grid 1,1 -dm_landau_Ez 0 -petscspace_degree 3 -petscspace_poly_tensor 1 -dm_landau_type p4est -info :tsadapt -dm_landau_ion_masses 2 -dm_landau_ion_charges 1 -dm_landau_thermal_temps 5,5 -dm_landau_n 2,2 -dm_landau_n_0 5e19 -ts_monitor -snes_rtol 1.e-10 -snes_stol 1.e-14 -snes_monitor -snes_converged_reason -snes_max_it 10 -ts_type arkimex -ts_arkimex_type 1bee -ts_max_snes_failures -1 -ts_rtol 1e-3 -ts_dt 1.e-1 -ts_max_time 1 -ts_adapt_clip .5,1.25 -ts_max_steps 2 -ts_adapt_scale_solve_failed 0.75 -ts_adapt_time_step_increase_delay 5 -pc_type lu -ksp_type preonly -dm_landau_amr_levels_max 3,2 -ex2_impurity_source_type pulse -ex2_pulse_start_time 1e-1 -ex2_pulse_width_time 10 -ex2_pulse_rate 1e-2 -ex2_t_cold .05 -ex2_plot_dt 1e-1 -dm_refine 1 -dm_landau_device_type cuda -dm_landau_gpu_assembly true -dm_mat_type aijcusparse -dm_vec_type cuda
771: TEST*/