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