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