Lines Matching refs:ctx
88 LandauCtx *ctx = (LandauCtx *)a_ctx; in LandauFormJacobian_Internal() local
102 PetscAssertPointer(ctx, 5); in LandauFormJacobian_Internal()
104 PetscCheck(ctx->plex[0] != NULL, ctx->comm, PETSC_ERR_ARG_WRONG, "Plex not created"); in LandauFormJacobian_Internal()
105 PetscCall(PetscLogEventBegin(ctx->events[10], 0, 0, 0, 0)); in LandauFormJacobian_Internal()
106 PetscCall(DMGetDS(ctx->plex[0], &prob)); // same DS for all grids in LandauFormJacobian_Internal()
109 PetscCheck(ctx->gpu_assembly, ctx->comm, PETSC_ERR_ARG_WRONG, "maps but no GPU assembly"); in LandauFormJacobian_Internal()
111 PetscCheck(maps, ctx->comm, PETSC_ERR_ARG_WRONG, "empty GPU matrix container"); in LandauFormJacobian_Internal()
112 for (PetscInt i = 0; i < ctx->num_grids * ctx->batch_sz; i++) subJ[i] = NULL; in LandauFormJacobian_Internal()
114 PetscCheck(!ctx->gpu_assembly, ctx->comm, PETSC_ERR_ARG_WRONG, "No maps but GPU assembly"); in LandauFormJacobian_Internal()
115 for (PetscInt tid = 0; tid < ctx->batch_sz; tid++) { in LandauFormJacobian_Internal()
116 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) PetscCall(DMCreateMatrix(ctx->plex[grid], &… in LandauFormJacobian_Internal()
121 PetscCall(PetscFEGetQuadrature(ctx->fe[0], &quad)); in LandauFormJacobian_Internal()
123 PetscCall(PetscFEGetDimension(ctx->fe[0], &Nb)); in LandauFormJacobian_Internal()
124 …PetscCheck(Nq <= LANDAU_MAX_NQND, ctx->comm, PETSC_ERR_ARG_WRONG, "Order too high. Nq = %" PetscIn… in LandauFormJacobian_Internal()
125 …PetscCheck(Nb <= LANDAU_MAX_NQND, ctx->comm, PETSC_ERR_ARG_WRONG, "Order too high. Nb = %" PetscIn… in LandauFormJacobian_Internal()
127 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauFormJacobian_Internal()
129 PetscCheck(ctx->plex[grid] != NULL, ctx->comm, PETSC_ERR_ARG_WRONG, "Plex not created"); in LandauFormJacobian_Internal()
130 PetscCall(DMPlexGetHeightStratum(ctx->plex[grid], 0, &cStart, &cEnd)); in LandauFormJacobian_Internal()
133 PetscCall(PetscLogEventEnd(ctx->events[10], 0, 0, 0, 0)); in LandauFormJacobian_Internal()
138 PetscCall(PetscLogEventBegin(ctx->events[1], 0, 0, 0, 0)); in LandauFormJacobian_Internal()
139 for (PetscInt fieldA = 0; fieldA < ctx->num_species; fieldA++) { in LandauFormJacobian_Internal()
140 …Eq_m[fieldA] = ctx->Ez * ctx->t_0 * ctx->charges[fieldA] / (ctx->v_0 * ctx->masses[fieldA]); /* no… in LandauFormJacobian_Internal()
143 if (!ctx->gpu_assembly) { in LandauFormJacobian_Internal()
148 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauFormJacobian_Internal()
149 PetscCall(DMGetLocalSection(ctx->plex[grid], §ion[grid])); in LandauFormJacobian_Internal()
150 PetscCall(DMGetGlobalSection(ctx->plex[grid], &globsection[grid])); in LandauFormJacobian_Internal()
155 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) cellClosure_sz += Nb * Nf[grid] * numCells[… in LandauFormJacobian_Internal()
156 PetscCall(PetscMalloc1(cellClosure_sz * ctx->batch_sz, &cellClosure)); in LandauFormJacobian_Internal()
162 for (PetscInt b_id = 0; b_id < ctx->batch_sz; b_id++) { // OpenMP (once) in LandauFormJacobian_Internal()
163 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauFormJacobian_Internal()
167 PetscCall(DMGlobalToLocalBegin(ctx->plex[grid], globX, INSERT_VALUES, locX2)); in LandauFormJacobian_Internal()
168 PetscCall(DMGlobalToLocalEnd(ctx->plex[grid], globX, INSERT_VALUES, locX2)); in LandauFormJacobian_Internal()
169 PetscCall(DMPlexGetHeightStratum(ctx->plex[grid], 0, &cStart, &cEnd)); in LandauFormJacobian_Internal()
172 PetscCall(DMPlexVecGetClosure(ctx->plex[grid], section[grid], locX2, ei, NULL, &coef)); in LandauFormJacobian_Internal()
174 … PetscCall(DMPlexVecRestoreClosure(ctx->plex[grid], section[grid], locX2, ei, NULL, &coef)); in LandauFormJacobian_Internal()
180 …ctx->batch_sz, PETSC_COMM_SELF, PETSC_ERR_PLIB, "iteration wrong %" PetscCount_FMT " != cellClosur… in LandauFormJacobian_Internal()
188 if (ctx->jacobian_field_major_order) { // get data in batch ordering in LandauFormJacobian_Internal()
189 … PetscCall(VecScatterBegin(ctx->plex_batch, a_X, ctx->work_vec, INSERT_VALUES, SCATTER_FORWARD)); in LandauFormJacobian_Internal()
190 … PetscCall(VecScatterEnd(ctx->plex_batch, a_X, ctx->work_vec, INSERT_VALUES, SCATTER_FORWARD)); in LandauFormJacobian_Internal()
191 PetscCall(VecGetArrayReadAndMemType(ctx->work_vec, &xdata, &mtype)); in LandauFormJacobian_Internal()
195 …PetscCheck(mtype == PETSC_MEMTYPE_HOST || ctx->deviceType != LANDAU_CPU, ctx->comm, PETSC_ERR_ARG_… in LandauFormJacobian_Internal()
198 PetscCall(PetscLogEventEnd(ctx->events[1], 0, 0, 0, 0)); in LandauFormJacobian_Internal()
202 if (ctx->deviceType == LANDAU_KOKKOS) { in LandauFormJacobian_Internal()
204 …obian(ctx->plex, Nq, Nb, ctx->batch_sz, ctx->num_grids, numCells, Eq_m, cellClosure, xdata, &ctx->… in LandauFormJacobian_Internal()
206 SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONG, "-landau_device_type %s not built", "kokkos"); in LandauFormJacobian_Internal()
210 …1], elem_offset[LANDAU_MAX_GRIDS + 1], IPf_sz_glb, IPf_sz_tot, num_grids = ctx->num_grids, Nf[LAND… in LandauFormJacobian_Internal()
211 … *)ctx->SData_d.invJ, *xx = (PetscReal *)ctx->SData_d.x, *yy = (PetscReal *)ctx->SData_d.y, *zz = … in LandauFormJacobian_Internal()
212 … *nu_alpha = (PetscReal *)ctx->SData_d.alpha, *nu_beta = (PetscReal *)ctx->SData_d.beta, *inv… in LandauFormJacobian_Internal()
213 …GRIDS][LANDAU_MAX_GRIDS] = (PetscReal (*)[LANDAU_MAX_GRIDS][LANDAU_MAX_GRIDS])ctx->SData_d.lambdas; in LandauFormJacobian_Internal()
216 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauFormJacobian_Internal()
217 PetscCall(DMGetLocalSection(ctx->plex[grid], §ion[grid])); in LandauFormJacobian_Internal()
218 PetscCall(DMGetGlobalSection(ctx->plex[grid], &globsection[grid])); in LandauFormJacobian_Internal()
226 PetscInt nfloc = ctx->species_offset[grid + 1] - ctx->species_offset[grid]; in LandauFormJacobian_Internal()
232 IPf_sz_tot = IPf_sz_glb * ctx->batch_sz; in LandauFormJacobian_Internal()
234 PetscCall(PetscMalloc1(ctx->SData_d.coo_size, &coo_vals)); // allocate every time? in LandauFormJacobian_Internal()
240 PetscCall(PetscLogEventBegin(ctx->events[8], 0, 0, 0, 0)); in LandauFormJacobian_Internal()
243 …for (PetscInt tid = 0; tid < ctx->batch_sz * elem_offset[num_grids]; tid++) { … in LandauFormJacobian_Internal()
249 …const PetscInt loc_nip = numCells[grid] * Nq, loc_Nf = ctx->species_offset[grid + 1] - ctx->specie… in LandauFormJacobian_Internal()
250 …const PetscInt moffset = LAND_MOFFSET(b_id, grid, ctx->batch_sz, ctx->num_grids, ctx->mat_offset);… in LandauFormJacobian_Internal()
306 PetscCall(PetscLogEventEnd(ctx->events[8], 0, 0, 0, 0)); in LandauFormJacobian_Internal()
309 if (ctx->stage) ctx->times[LANDAU_F_DF] += (endtime - starttime); in LandauFormJacobian_Internal()
313 … for (PetscInt tid = 0; tid < ctx->batch_sz * elem_offset[num_grids]; tid++) { // for each element in LandauFormJacobian_Internal()
323 …const PetscInt loc_Nf = ctx->species_offset[grid + 1] - ctx->species_offset[grid], loc_elem = gl… in LandauFormJacobian_Internal()
324 …const PetscInt moffset = LAND_MOFFSET(b_id, grid, ctx->batch_sz, ctx->num_grids, ctx->mat_offset… in LandauFormJacobian_Internal()
330 PetscCall(PetscLogEventBegin(ctx->events[4], 0, 0, 0, 0)); in LandauFormJacobian_Internal()
332 PetscCall(PetscLogEventBegin(ctx->events[16], 0, 0, 0, 0)); in LandauFormJacobian_Internal()
349 …for (PetscInt grid_r = 0, f_off = 0, ipidx = 0; grid_r < ctx->num_grids; grid_r++, f_off = ctx->sp… in LandauFormJacobian_Internal()
360 if (ctx->use_relativistic_corrections) { in LandauFormJacobian_Internal()
361 LandauTensor3DRelativistic(vj, x, y, z, U, mask, C_0(ctx->v_0)); in LandauFormJacobian_Internal()
407 … for (PetscInt fieldA = 0, f_off = ctx->species_offset[grid]; fieldA < loc_Nf; ++fieldA) { in LandauFormJacobian_Internal()
414 …for (PetscInt fieldA = 0, f_off = ctx->species_offset[grid]; fieldA < loc_Nf; ++fieldA) gg2[fieldA… in LandauFormJacobian_Internal()
466 PetscCall(PetscLogEventEnd(ctx->events[4], 0, 0, 0, 0)); in LandauFormJacobian_Internal()
468 PetscCall(PetscLogEventEnd(ctx->events[16], 0, 0, 0, 0)); in LandauFormJacobian_Internal()
472 if (ctx->stage) ctx->times[LANDAU_KERNEL] += (endtime - starttime); in LandauFormJacobian_Internal()
477 PetscCall(PetscLogEventBegin(ctx->events[6], 0, 0, 0, 0)); in LandauFormJacobian_Internal()
478 PetscCall(DMPlexGetHeightStratum(ctx->plex[grid], 0, &cStart, NULL)); in LandauFormJacobian_Internal()
479 …PetscCall(DMPlexMatSetClosure(ctx->plex[grid], section[grid], globsection[grid], subJ[LAND_PACK_ID… in LandauFormJacobian_Internal()
480 PetscCall(PetscLogEventEnd(ctx->events[6], 0, 0, 0, 0)); in LandauFormJacobian_Internal()
484 …ctx->SData_d.coo_elem_offsets, *coo_elem_fullNb = (LandauIdx *)ctx->SData_d.coo_elem_fullNb, (*coo… in LandauFormJacobian_Internal()
539 PetscCall(PetscPrintf(ctx->comm, "CPU Element matrix\n")); in LandauFormJacobian_Internal()
541 …for (PetscInt f = 0; f < totDim; ++f) PetscCall(PetscPrintf(ctx->comm, " %12.5e", (double)PetscRea… in LandauFormJacobian_Internal()
542 PetscCall(PetscPrintf(ctx->comm, "\n")); in LandauFormJacobian_Internal()
553 for (PetscInt b_id = 0; b_id < ctx->batch_sz; b_id++) { // OpenMP in LandauFormJacobian_Internal()
554 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauFormJacobian_Internal()
555 …const PetscInt moffset = LAND_MOFFSET(b_id, grid, ctx->batch_sz, ctx->num_grids, ctx->mat_offs… in LandauFormJacobian_Internal()
589 …CreateVMeshes(MPI_Comm comm_self, const PetscInt dim, const char prefix[], LandauCtx *ctx, DM pack) in LandauDMCreateVMeshes() argument
594 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauDMCreateVMeshes()
595 PetscReal par_radius = ctx->radius_par[grid], perp_radius = ctx->radius_perp[grid]; in LandauDMCreateVMeshes()
596 if (!ctx->sphere && !ctx->simplex) { // 2 or 3D (only 3D option) in LandauDMCreateVMeshes()
604 …xCreateBoxMesh(comm_self, dim, PETSC_FALSE, ctx->cells0, lo, hi, periodicity, PETSC_TRUE, 0, PETSC… in LandauDMCreateVMeshes()
605 …PetscCall(DMLocalizeCoordinates(ctx->plex[grid])); … in LandauDMCreateVMeshes()
606 if (dim == 3) PetscCall(PetscObjectSetName((PetscObject)ctx->plex[grid], "cube")); in LandauDMCreateVMeshes()
607 else PetscCall(PetscObjectSetName((PetscObject)ctx->plex[grid], "half-plane")); in LandauDMCreateVMeshes()
610 PetscCall(PetscStrlen(ctx->filename, &len)); in LandauDMCreateVMeshes()
617 …PetscCall(DMPlexCreateFromFile(comm_self, ctx->filename, "plexland.c", PETSC_TRUE, &ctx->plex[grid… in LandauDMCreateVMeshes()
618 PetscCall(DMPlexOrient(ctx->plex[grid])); in LandauDMCreateVMeshes()
619 PetscCall(DMGetCoordinatesLocal(ctx->plex[grid], &coords)); in LandauDMCreateVMeshes()
624 x[i + 0] *= ctx->radius_perp[grid]; in LandauDMCreateVMeshes()
625 x[i + 1] *= ctx->radius_par[grid]; in LandauDMCreateVMeshes()
628 PetscCall(PetscObjectSetName((PetscObject)ctx->plex[grid], ctx->filename)); in LandauDMCreateVMeshes()
629 …PetscCall(PetscInfo(ctx->plex[grid], "%" PetscInt_FMT ") Read %s mesh file (%s)\n", grid, ctx->fil… in LandauDMCreateVMeshes()
630 PetscCall(DMViewFromOptions(ctx->plex[grid], NULL, str)); in LandauDMCreateVMeshes()
632 PetscInt numCells = ctx->simplex ? 12 : 6, cell_size = ctx->simplex ? 3 : 4, j; in LandauDMCreateVMeshes()
656 const PetscInt *pcell = (const PetscInt *)(ctx->simplex ? &cellsS[0][0] : &cellsT[0][0]); in LandauDMCreateVMeshes()
658 PetscReal rad = ctx->radius[grid]; in LandauDMCreateVMeshes()
667 coords[j++][1] = -rad * ctx->sphere_inner_radius_90degree[grid]; in LandauDMCreateVMeshes()
668 coords[j][0] = rad * ctx->sphere_inner_radius_45degree[grid] * 0.707106781186548; in LandauDMCreateVMeshes()
669 coords[j++][1] = -rad * ctx->sphere_inner_radius_45degree[grid] * 0.707106781186548; in LandauDMCreateVMeshes()
670 coords[j][0] = rad * ctx->sphere_inner_radius_90degree[grid]; in LandauDMCreateVMeshes()
672 coords[j][0] = rad * ctx->sphere_inner_radius_45degree[grid] * 0.707106781186548; in LandauDMCreateVMeshes()
673 coords[j++][1] = rad * ctx->sphere_inner_radius_45degree[grid] * 0.707106781186548; in LandauDMCreateVMeshes()
675 coords[j++][1] = rad * ctx->sphere_inner_radius_90degree[grid]; in LandauDMCreateVMeshes()
678 …istPetsc(comm_self, 2, numCells, numVerts, cell_size, ctx->interpolate, pcell, 2, flatCoords, &ctx… in LandauDMCreateVMeshes()
679 PetscCall(PetscObjectSetName((PetscObject)ctx->plex[grid], "semi-circle")); in LandauDMCreateVMeshes()
680 …PetscCall(PetscInfo(ctx->plex[grid], "\t%" PetscInt_FMT ") Make circle %s mesh\n", grid, ctx->simp… in LandauDMCreateVMeshes()
683 …PetscCheck(dim == 3 && ctx->sphere && !ctx->simplex, ctx->comm, PETSC_ERR_ARG_WRONG, "not: dim == … in LandauDMCreateVMeshes()
684 …PetscReal rad = ctx->radius[grid], inner_rad = rad * ctx->sphere_inner_radius_90degree[grid],… in LandauDMCreateVMeshes()
716 …lf, 3, numCells, numVerts, cell_size, ctx->interpolate, (const PetscInt *)cells, 3, (const PetscRe… in LandauDMCreateVMeshes()
717 PetscCall(PetscObjectSetName((PetscObject)ctx->plex[grid], "cubed sphere")); in LandauDMCreateVMeshes()
718 …PetscCall(PetscInfo(ctx->plex[grid], "\t%" PetscInt_FMT ") Make cubed sphere %s mesh\n", grid, ctx… in LandauDMCreateVMeshes()
720 PetscCall(DMSetOptionsPrefix(ctx->plex[grid], prefix)); in LandauDMCreateVMeshes()
721 PetscCall(DMSetFromOptions(ctx->plex[grid])); in LandauDMCreateVMeshes()
728 PetscOptionsBegin(ctx->comm, prefix, "Mesh conversion options", "DMPLEX"); in LandauDMCreateVMeshes()
732 ctx->use_p4est = PETSC_TRUE; /* flag for Forest */ in LandauDMCreateVMeshes()
733 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauDMCreateVMeshes()
737 PetscCall(DMConvert(ctx->plex[grid], convType, &dmforest)); in LandauDMCreateVMeshes()
738 PetscCheck(dmforest, ctx->comm, PETSC_ERR_PLIB, "Convert failed?"); in LandauDMCreateVMeshes()
741 PetscCheck(isForest, ctx->comm, PETSC_ERR_PLIB, "Converted to non Forest?"); in LandauDMCreateVMeshes()
742 PetscCall(DMDestroy(&ctx->plex[grid])); in LandauDMCreateVMeshes()
743 ctx->plex[grid] = dmforest; // Forest for adaptivity in LandauDMCreateVMeshes()
745 } else ctx->use_p4est = PETSC_FALSE; /* flag for Forest */ in LandauDMCreateVMeshes()
749 PetscCall(DMSetApplicationContext(pack, ctx)); in LandauDMCreateVMeshes()
753 …c PetscErrorCode SetupDS(DM pack, PetscInt dim, PetscInt grid, const char prefix[], LandauCtx *ctx) in SetupDS() argument
760 for (ii = ctx->species_offset[grid], i0 = 0; ii < ctx->species_offset[grid + 1]; ii++, i0++) { in SetupDS()
764 …PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, ctx->simplex, prefix, PETSC_DECIDE, &ctx->… in SetupDS()
765 PetscCall(PetscObjectSetName((PetscObject)ctx->fe[ii], buf)); in SetupDS()
766 PetscCall(DMSetField(ctx->plex[grid], i0, NULL, (PetscObject)ctx->fe[ii])); in SetupDS()
768 PetscCall(DMCreateDS(ctx->plex[grid])); in SetupDS()
769 PetscCall(DMGetLocalSection(ctx->plex[grid], §ion)); in SetupDS()
770 …for (PetscInt ii = ctx->species_offset[grid], i0 = 0; ii < ctx->species_offset[grid + 1]; ii++, i0… in SetupDS()
851 LandauCtx *ctx = (LandauCtx *)actx; in DMPlexLandauAddMaxwellians() local
858 if (!ctx) PetscCall(DMGetApplicationContext(dm, &ctx)); in DMPlexLandauAddMaxwellians()
859 …for (PetscInt ii = ctx->species_offset[grid], i0 = 0; ii < ctx->species_offset[grid + 1]; ii++, i0… in DMPlexLandauAddMaxwellians()
861 data[i0].v_0 = ctx->v_0; // v_0 same for all grids in DMPlexLandauAddMaxwellians()
862 data[i0].kT_m = ctx->k * temps[ii] / ctx->masses[ii]; /* kT/m */ in DMPlexLandauAddMaxwellians()
867 data[0].shift = ctx->electronShift; in DMPlexLandauAddMaxwellians()
894 LandauCtx *ctx = (LandauCtx *)actx; in LandauSetInitialCondition() local
897 if (!ctx) PetscCall(DMGetApplicationContext(dm, &ctx)); in LandauSetInitialCondition()
899 …PetscCall(DMPlexLandauAddMaxwellians(dm, X, 0.0, ctx->thermal_temps, ctx->n, grid, b_id, n_batch, … in LandauSetInitialCondition()
907 …daptToleranceFEM(PetscFE fem, Vec sol, PetscInt type, PetscInt grid, LandauCtx *ctx, DM *newForest) in adaptToleranceFEM() argument
917 forest = ctx->plex[grid]; in adaptToleranceFEM()
922 PetscCheck(isForest, ctx->comm, PETSC_ERR_ARG_WRONG, "! Forest"); in adaptToleranceFEM()
928 …PetscCheck(Nq <= LANDAU_MAX_NQND, ctx->comm, PETSC_ERR_ARG_WRONG, "Order too high. Nq = %" PetscIn… in adaptToleranceFEM()
929 PetscCall(PetscFEGetDimension(ctx->fe[0], &Nb)); in adaptToleranceFEM()
931 …PetscCheck(Nb2[0] == Nb, ctx->comm, PETSC_ERR_ARG_WRONG, " Nb = %" PetscInt_FMT " != Nb (%" PetscI… in adaptToleranceFEM()
932 …PetscCheck(Nb <= LANDAU_MAX_NQND, ctx->comm, PETSC_ERR_ARG_WRONG, "Order too high. Nb = %" PetscIn… in adaptToleranceFEM()
981 …if (ctx->re_radius > PETSC_SQRT_MACHINE_EPSILON && (z < -PETSC_MACHINE_EPSILON * 10. || z > ctx->r… in adaptToleranceFEM()
982 } else if (type == 1 && (z > ctx->vperp0_radius1 || z < -ctx->vperp0_radius1)) { in adaptToleranceFEM()
985 } else if (type == 3 && (z > ctx->vperp0_radius2 || z < -ctx->vperp0_radius2)) { in adaptToleranceFEM()
989 …if (x < PETSC_MACHINE_EPSILON * 10. && (type != 0 || ctx->re_radius > PETSC_SQRT_MACHINE_EPSILON))… in adaptToleranceFEM()
1014 static PetscErrorCode adapt(PetscInt grid, LandauCtx *ctx, Vec *uu) in adapt() argument
1019 …= {(grid == 0) ? ctx->numRERefine : 0, (grid == 0) ? ctx->nZRefine1 : 0, ctx->numAMRRefine[grid], … in adapt()
1023 PetscCall(adaptToleranceFEM(ctx->fe[0], *uu, type, grid, ctx, &newForest)); in adapt()
1025 PetscCall(DMDestroy(&ctx->plex[grid])); in adapt()
1028 PetscCall(LandauSetInitialCondition(newForest, *uu, grid, 0, 1, ctx)); in adapt()
1029 ctx->plex[grid] = newForest; in adapt()
1040 static PetscErrorCode makeLambdas(LandauCtx *ctx) in makeLambdas() argument
1043 for (PetscInt gridi = 0; gridi < ctx->num_grids; gridi++) { in makeLambdas()
1044 PetscInt iii = ctx->species_offset[gridi]; in makeLambdas()
1045 PetscReal Ti_ev = (ctx->thermal_temps[iii] / 1.1604525e7) * 1000; // convert (back) to eV in makeLambdas()
1046 PetscReal ni = ctx->n[iii] * ctx->n_0; in makeLambdas()
1047 for (PetscInt gridj = gridi; gridj < ctx->num_grids; gridj++) { in makeLambdas()
1048 PetscInt jjj = ctx->species_offset[gridj]; in makeLambdas()
1049 PetscReal Zj = ctx->charges[jjj] / 1.6022e-19; in makeLambdas()
1052 …ctx->lambdas[gridi][gridj] = 23.5 - PetscLogReal(PetscSqrtReal(ni) * PetscPowReal(Ti_ev, -1.25)) -… in makeLambdas()
1055 …ctx->lambdas[gridi][gridj] = ctx->lambdas[gridj][gridi] = 23 - PetscLogReal(PetscSqrtReal(ni) * Zj… in makeLambdas()
1057 …ctx->lambdas[gridi][gridj] = ctx->lambdas[gridj][gridi] = 24 - PetscLogReal(PetscSqrtReal(ni) / Ti… in makeLambdas()
1061 PetscReal mui = ctx->masses[iii] / 1.6720e-27, Zi = ctx->charges[iii] / 1.6022e-19; in makeLambdas()
1062 …PetscReal Tj_ev = (ctx->thermal_temps[jjj] / 1.1604525e7) * 1000; // convert (back) to … in makeLambdas()
1063 PetscReal muj = ctx->masses[jjj] / 1.6720e-27; in makeLambdas()
1064 PetscReal nj = ctx->n[jjj] * ctx->n_0; in makeLambdas()
1065 …ctx->lambdas[gridi][gridj] = ctx->lambdas[gridj][gridi] = 23 - PetscLogReal(Zi * Zj * (mui + muj) … in makeLambdas()
1073 static PetscErrorCode ProcessOptions(LandauCtx *ctx, const char prefix[]) in ProcessOptions() argument
1081 PetscCall(DMCreate(ctx->comm, &dummy)); in ProcessOptions()
1083 ctx->verbose = 1; // should be 0 for silent compliance in ProcessOptions()
1084 ctx->batch_sz = 1; in ProcessOptions()
1085 ctx->batch_view_idx = 0; in ProcessOptions()
1086 ctx->interpolate = PETSC_TRUE; in ProcessOptions()
1087 ctx->gpu_assembly = PETSC_TRUE; in ProcessOptions()
1088 ctx->norm_state = 0; in ProcessOptions()
1089 ctx->electronShift = 0; in ProcessOptions()
1090 ctx->M = NULL; in ProcessOptions()
1091 ctx->J = NULL; in ProcessOptions()
1093 ctx->sphere = PETSC_FALSE; in ProcessOptions()
1094 ctx->map_sphere = PETSC_TRUE; in ProcessOptions()
1095 ctx->use_p4est = PETSC_FALSE; in ProcessOptions()
1096 ctx->simplex = PETSC_FALSE; in ProcessOptions()
1098 ctx->radius[grid] = 5.; /* thermal radius (velocity) */ in ProcessOptions()
1099 ctx->radius_perp[grid] = 5.; /* thermal radius (velocity) */ in ProcessOptions()
1100 ctx->radius_par[grid] = 5.; /* thermal radius (velocity) */ in ProcessOptions()
1101 ctx->numAMRRefine[grid] = 0; in ProcessOptions()
1102 ctx->postAMRRefine[grid] = 0; in ProcessOptions()
1103 ctx->species_offset[grid + 1] = 1; // one species default in ProcessOptions()
1105 ctx->plex[grid] = NULL; /* cache as expensive to Convert */ in ProcessOptions()
1107 ctx->species_offset[0] = 0; in ProcessOptions()
1108 ctx->re_radius = 0.; in ProcessOptions()
1109 ctx->vperp0_radius1 = 0; in ProcessOptions()
1110 ctx->vperp0_radius2 = 0; in ProcessOptions()
1111 ctx->nZRefine1 = 0; in ProcessOptions()
1112 ctx->nZRefine2 = 0; in ProcessOptions()
1113 ctx->numRERefine = 0; in ProcessOptions()
1116 ctx->charges[0] = -1; /* electron charge (MKS) */ in ProcessOptions()
1117 ctx->masses[0] = 1 / 1835.469965278441013; /* temporary value in proton mass */ in ProcessOptions()
1118 ctx->n[0] = 1; in ProcessOptions()
1119 ctx->v_0 = 1; /* thermal velocity, we could start with a scale != 1 */ in ProcessOptions()
1120 ctx->thermal_temps[0] = 1; in ProcessOptions()
1122 ctx->epsilon0 = 8.8542e-12; /* permittivity of free space (MKS) F/m */ in ProcessOptions()
1123 ctx->k = 1.38064852e-23; /* Boltzmann constant (MKS) J/K */ in ProcessOptions()
1124 ctx->n_0 = 1.e20; /* typical plasma n, but could set it to 1 */ in ProcessOptions()
1125 ctx->Ez = 0; in ProcessOptions()
1126 for (PetscInt grid = 0; grid < LANDAU_NUM_TIMERS; grid++) ctx->times[grid] = 0; in ProcessOptions()
1127 for (PetscInt ii = 0; ii < LANDAU_DIM; ii++) ctx->cells0[ii] = 2; in ProcessOptions()
1128 if (LANDAU_DIM == 2) ctx->cells0[0] = 1; in ProcessOptions()
1129 ctx->use_matrix_mass = PETSC_FALSE; in ProcessOptions()
1130 ctx->use_relativistic_corrections = PETSC_FALSE; in ProcessOptions()
1131 …ctx->use_energy_tensor_trick = PETSC_FALSE; /* Use Eero's trick for energy conservation v -… in ProcessOptions()
1132 ctx->SData_d.w = NULL; in ProcessOptions()
1133 ctx->SData_d.x = NULL; in ProcessOptions()
1134 ctx->SData_d.y = NULL; in ProcessOptions()
1135 ctx->SData_d.z = NULL; in ProcessOptions()
1136 ctx->SData_d.invJ = NULL; in ProcessOptions()
1137 ctx->jacobian_field_major_order = PETSC_FALSE; in ProcessOptions()
1138 ctx->SData_d.coo_elem_offsets = NULL; in ProcessOptions()
1139 ctx->SData_d.coo_elem_point_offsets = NULL; in ProcessOptions()
1140 ctx->SData_d.coo_elem_fullNb = NULL; in ProcessOptions()
1141 ctx->SData_d.coo_size = 0; in ProcessOptions()
1142 …PetscOptionsBegin(ctx->comm, prefix, "Options for Fokker-Plank-Landau collision operator", "none"); in ProcessOptions()
1144 ctx->deviceType = LANDAU_KOKKOS; in ProcessOptions()
1145 PetscCall(PetscStrncpy(ctx->filename, "kokkos", sizeof(ctx->filename))); in ProcessOptions()
1147 ctx->deviceType = LANDAU_CPU; in ProcessOptions()
1148 PetscCall(PetscStrncpy(ctx->filename, "cpu", sizeof(ctx->filename))); in ProcessOptions()
1150 …e_type", "Use kernels on 'cpu' 'kokkos'", "plexland.c", ctx->filename, ctx->filename, sizeof(ctx->… in ProcessOptions()
1151 PetscCall(PetscStrcmp("cpu", ctx->filename, &flg)); in ProcessOptions()
1153 ctx->deviceType = LANDAU_CPU; in ProcessOptions()
1155 PetscCall(PetscStrcmp("kokkos", ctx->filename, &flg)); in ProcessOptions()
1156 PetscCheck(flg, ctx->comm, PETSC_ERR_ARG_WRONG, "-dm_landau_device_type %s", ctx->filename); in ProcessOptions()
1157 ctx->deviceType = LANDAU_KOKKOS; in ProcessOptions()
1159 ctx->filename[0] = '\0'; in ProcessOptions()
1160 …ndau_filename", "file to read mesh from", "plexland.c", ctx->filename, ctx->filename, sizeof(ctx->… in ProcessOptions()
1161 …ctron_shift", "Shift in thermal velocity of electrons", "none", ctx->electronShift, &ctx->electron… in ProcessOptions()
1162 …t("-dm_landau_verbose", "Level of verbosity output", "plexland.c", ctx->verbose, &ctx->verbose, NU… in ProcessOptions()
1163 …"-dm_landau_batch_size", "Number of 'vertices' to batch", "ex2.c", ctx->batch_sz, &ctx->batch_sz, … in ProcessOptions()
1164 …BATCH_SZ >= ctx->batch_sz, ctx->comm, PETSC_ERR_ARG_WRONG, "LANDAU_MAX_BATCH_SZ %d < ctx->batch_sz… in ProcessOptions()
1165 …_idx", "Index of batch for diagnostics like plotting", "ex2.c", ctx->batch_view_idx, &ctx->batch_v… in ProcessOptions()
1166 …ctx->batch_view_idx < ctx->batch_sz, ctx->comm, PETSC_ERR_ARG_WRONG, "-ctx->batch_view_idx %" Pets… in ProcessOptions()
1167 …electric field in unites of Conner-Hastie critical field", "plexland.c", ctx->Ez, &ctx->Ez, NULL)); in ProcessOptions()
1168 …andau_n_0", "Normalization constant for number density", "plexland.c", ctx->n_0, &ctx->n_0, NULL)); in ProcessOptions()
1169 … but slightly fragile MATAXPY to add mass term", "plexland.c", ctx->use_matrix_mass, &ctx->use_mat… in ProcessOptions()
1170 …rections", "Use relativistic corrections", "plexland.c", ctx->use_relativistic_corrections, &ctx->… in ProcessOptions()
1171 …nsBool("-dm_landau_simplex", "Use simplex elements", "plexland.c", ctx->simplex, &ctx->simplex, NU… in ProcessOptions()
1172 …"use sphere/semi-circle domain instead of rectangle", "plexland.c", ctx->sphere, &ctx->sphere, NUL… in ProcessOptions()
1173 …to sphere/semi-circle domain instead of rectangle", "plexland.c", ctx->map_sphere, &ctx->map_spher… in ProcessOptions()
1174 …if (LANDAU_DIM == 2 && ctx->use_relativistic_corrections) ctx->use_relativistic_corrections = PETS… in ProcessOptions()
1175 …serve energy with relativistic corrections and Q1 elements", "plexland.c", ctx->use_energy_tensor_… in ProcessOptions()
1176 &ctx->use_energy_tensor_trick, NULL)); in ProcessOptions()
1180 ctx->thermal_temps[ii] = 1; in ProcessOptions()
1181 ctx->charges[ii] = 1; in ProcessOptions()
1182 ctx->masses[ii] = 1; in ProcessOptions()
1183 ctx->n[ii] = 1; in ProcessOptions()
1186 …i_0,i_1,...] in keV (must be set to set number of species)", "plexland.c", ctx->thermal_temps, &nt… in ProcessOptions()
1187 …PetscCheck(flg, ctx->comm, PETSC_ERR_ARG_WRONG, "-dm_landau_thermal_temps ,t1,t2,.. must be provid… in ProcessOptions()
1189 ctx->num_species = nt; in ProcessOptions()
1190 …for (ii = 0; ii < ctx->num_species; ii++) ctx->thermal_temps[ii] *= 1.1604525e7; /* convert to Kel… in ProcessOptions()
1192 …of each species in units of proton mass [i_0=2,i_1=40...]", "plexland.c", &ctx->masses[1], &nm, &f… in ProcessOptions()
1193 …flg || nm == ctx->num_species - 1, ctx->comm, PETSC_ERR_ARG_WRONG, "num ion masses %" PetscInt_FMT… in ProcessOptions()
1195 …ay("-dm_landau_n", "Number density of each species = n_s * n_0", "plexland.c", ctx->n, &nm, &flg)); in ProcessOptions()
1196 …k(!flg || nm == ctx->num_species, ctx->comm, PETSC_ERR_ARG_WRONG, "wrong num n: %" PetscInt_FMT " … in ProcessOptions()
1197 …for (ii = 0; ii < LANDAU_MAX_SPECIES; ii++) ctx->masses[ii] *= 1.6720e-27; /* scale by proton mass… in ProcessOptions()
1198 …ctx->masses[0] = 9.10938356e-31; /* electron mass kg (sh… in ProcessOptions()
1200 …each species in units of proton charge [i_0=2,i_1=18,...]", "plexland.c", &ctx->charges[1], &nc, &… in ProcessOptions()
1201 …scCheck(nc == ctx->num_species - 1, ctx->comm, PETSC_ERR_ARG_WRONG, "num charges %" PetscInt_FMT "… in ProcessOptions()
1202 …for (ii = 0; ii < LANDAU_MAX_SPECIES; ii++) ctx->charges[ii] *= 1.6022e-19; /* electron/proton cha… in ProcessOptions()
1207 ctx->num_grids = nt; in ProcessOptions()
1208 for (ii = nt = 0; ii < ctx->num_grids; ii++) nt += num_species_grid[ii]; in ProcessOptions()
1209 …ctx->num_species == nt, ctx->comm, PETSC_ERR_ARG_WRONG, "-dm_landau_num_species_grid: sum %" Petsc… in ProcessOptions()
1210 ctx->num_grids, LANDAU_MAX_GRIDS); in ProcessOptions()
1212 if (ctx->num_species > LANDAU_MAX_GRIDS) { in ProcessOptions()
1214 num_species_grid[1] = ctx->num_species - 1; in ProcessOptions()
1215 ctx->num_grids = 2; in ProcessOptions()
1217 ctx->num_grids = ctx->num_species; in ProcessOptions()
1218 for (ii = 0; ii < ctx->num_grids; ii++) num_species_grid[ii] = 1; in ProcessOptions()
1221 …for (ctx->species_offset[0] = ii = 0; ii < ctx->num_grids; ii++) ctx->species_offset[ii + 1] = ctx… in ProcessOptions()
1222 …ctx->species_offset[ctx->num_grids] == ctx->num_species, ctx->comm, PETSC_ERR_ARG_WRONG, "ctx->spe… in ProcessOptions()
1223 ctx->num_species); in ProcessOptions()
1224 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in ProcessOptions()
1225 …PetscInt iii = ctx->species_offset[grid]; // … in ProcessOptions()
1226 …ctx->thermal_speed[grid] = PetscSqrtReal(ctx->k * ctx->thermal_temps[iii] / ctx->masses[iii]); /* … in ProcessOptions()
1232 …for (PetscInt gridj = 0; gridj < LANDAU_MAX_GRIDS; gridj++) ctx->lambdas[gridj][grid] = lnLam; /* … in ProcessOptions()
1235 PetscCall(makeLambdas(ctx)); in ProcessOptions()
1240 …PetscCheck(non_dim_grid >= 0 && non_dim_grid < ctx->num_species, ctx->comm, PETSC_ERR_ARG_WRONG, "… in ProcessOptions()
1241 …ctx->v_0 = ctx->thermal_speed[non_dim_grid]; /* arbitrary units for non dimensionalization: global… in ProcessOptions()
1242 ctx->m_0 = ctx->masses[non_dim_grid]; /* arbitrary reference mass, electrons */ in ProcessOptions()
1243 …ctx->t_0 = 8 * PETSC_PI * PetscSqr(ctx->epsilon0 * ctx->m_0 / PetscSqr(ctx->charges[non_dim_grid])… in ProcessOptions()
1246 …s", "Phase space size in units of thermal velocity of grid", "plexland.c", ctx->radius, &nt, &flg)… in ProcessOptions()
1248 …>= ctx->num_grids, ctx->comm, PETSC_ERR_ARG_WRONG, "-dm_landau_domain_radius: given %" PetscInt_FM… in ProcessOptions()
1249 while (nt--) ctx->radius_par[nt] = ctx->radius_perp[nt] = ctx->radius[nt]; in ProcessOptions()
1252 … velocity domain size in units of thermal velocity of grid", "plexland.c", ctx->radius_par, &nt, &… in ProcessOptions()
1253 …>= ctx->num_grids, ctx->comm, PETSC_ERR_ARG_WRONG, "-dm_landau_domain_max_par: given %" PetscInt_F… in ProcessOptions()
1254 … velocity domain size in units of thermal velocity of grid", "plexland.c", ctx->radius_perp, &nt, … in ProcessOptions()
1255 …= ctx->num_grids, ctx->comm, PETSC_ERR_ARG_WRONG, "-dm_landau_domain_max_perp: given %" PetscInt_F… in ProcessOptions()
1257 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in ProcessOptions()
1258 …if (flg && ctx->radius[grid] <= 0) { /* negative is ratio of c - need to set par and perp with thi… in ProcessOptions()
1259 if (ctx->radius[grid] == 0) ctx->radius[grid] = 0.75; in ProcessOptions()
1260 else ctx->radius[grid] = -ctx->radius[grid]; in ProcessOptions()
1261 …ctx->radius[grid] = ctx->radius[grid] * SPEED_OF_LIGHT / ctx->v_0; // use any species on grid to n… in ProcessOptions()
1262 …(dummy, "Change domain radius to %g for grid %" PetscInt_FMT "\n", (double)ctx->radius[grid], grid… in ProcessOptions()
1264 …ctx->radius[grid] *= ctx->thermal_speed[grid] / ctx->v_0; // scale domain by thermal radius r… in ProcessOptions()
1265 …ctx->radius_perp[grid] *= ctx->thermal_speed[grid] / ctx->v_0; // scale domain by thermal radius r… in ProcessOptions()
1266 …ctx->radius_par[grid] *= ctx->thermal_speed[grid] / ctx->v_0; // scale domain by thermal radius r… in ProcessOptions()
1271 …f refinement around origin, after (RE) refinements along z", "plexland.c", ctx->numAMRRefine, &nt,… in ProcessOptions()
1272 …nt >= ctx->num_grids, ctx->comm, PETSC_ERR_ARG_WRONG, "-dm_landau_amr_levels_max: given %" PetscIn… in ProcessOptions()
1274 …t_refine", "Number of levels to uniformly refine after AMR", "plexland.c", ctx->postAMRRefine, &nt… in ProcessOptions()
1275 …for (ii = 1; ii < ctx->num_grids; ii++) ctx->postAMRRefine[ii] = ctx->postAMRRefine[0]; // all gri… in ProcessOptions()
1276 … "Number of levels to refine along v_perp=0, z>0", "plexland.c", ctx->numRERefine, &ctx->numRERefi… in ProcessOptions()
1277 …els to refine along v_perp=0 before origin refine", "plexland.c", ctx->nZRefine1, &ctx->nZRefine1,… in ProcessOptions()
1278 …vels to refine along v_perp=0 after origin refine", "plexland.c", ctx->nZRefine2, &ctx->nZRefine2,… in ProcessOptions()
1279 …to refine on positive (z>0) r=0 axis for runaways", "plexland.c", ctx->re_radius, &ctx->re_radius,… in ProcessOptions()
1280 …locity range to refine r=0 axis (for electrons)", "plexland.c", ctx->vperp0_radius1, &ctx->vperp0_… in ProcessOptions()
1281 …efine r=0 axis (for electrons) after origin AMR", "plexland.c", ctx->vperp0_radius2, &ctx->vperp0_… in ProcessOptions()
1283 if (ctx->sphere || ctx->simplex) { in ProcessOptions()
1284 ctx->sphere_uniform_normal = PETSC_FALSE; in ProcessOptions()
1285 …ticles per cell with Maxwellians (not used)", "plexland.c", ctx->sphere_uniform_normal, &ctx->sphe… in ProcessOptions()
1286 if (!ctx->sphere_uniform_normal) { // true in ProcessOptions()
1288 …le", "Scaling of radius for inner circle on 90 degree grid", "plexland.c", ctx->sphere_inner_radiu… in ProcessOptions()
1289 if (flg && nt < ctx->num_grids) { in ProcessOptions()
1290 …for (PetscInt grid = nt; grid < ctx->num_grids; grid++) ctx->sphere_inner_radius_90degree[grid] = … in ProcessOptions()
1292 if (ctx->sphere && !ctx->simplex && LANDAU_DIM == 3) { in ProcessOptions()
1293 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) ctx->sphere_inner_radius_90degree[grid] = 0… in ProcessOptions()
1296 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) ctx->sphere_inner_radius_90degree[grid] = 0… in ProcessOptions()
1298 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) ctx->sphere_inner_radius_90degree[grid] = 0… in ProcessOptions()
1303 …le", "Scaling of radius for inner circle on 45 degree grid", "plexland.c", ctx->sphere_inner_radiu… in ProcessOptions()
1304 if (flg && nt < ctx->num_grids) { in ProcessOptions()
1305 …for (PetscInt grid = nt; grid < ctx->num_grids; grid++) ctx->sphere_inner_radius_45degree[grid] = … in ProcessOptions()
1308 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) ctx->sphere_inner_radius_45degree[grid] = 0… in ProcessOptions()
1310 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) ctx->sphere_inner_radius_45degree[grid] = 0… in ProcessOptions()
1313 …ctx->sphere) PetscCall(PetscInfo(ctx->plex[0], "sphere : , 45 degree scaling = %g; 90 degree scali… in ProcessOptions()
1315 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in ProcessOptions()
1316 switch (ctx->numAMRRefine[grid]) { in ProcessOptions()
1323 ctx->sphere_inner_radius_90degree[grid] = 0.40; in ProcessOptions()
1324 ctx->sphere_inner_radius_45degree[grid] = 0.45; in ProcessOptions()
1326 ctx->sphere_inner_radius_45degree[grid] = 0.25; in ProcessOptions()
1333 …um_cells", "Number of cells in each dimension of base grid", "plexland.c", ctx->cells0, &nt, &flg)… in ProcessOptions()
1337 …_landau_gpu_assembly", "Assemble Jacobian on GPU", "plexland.c", ctx->gpu_assembly, &ctx->gpu_asse… in ProcessOptions()
1338 … or block diagonal, ordering (DEPRECATED)", "plexland.c", ctx->jacobian_field_major_order, &ctx->j… in ProcessOptions()
1339 …if (ctx->jacobian_field_major_order) PetscCheck(ctx->gpu_assembly, ctx->comm, PETSC_ERR_ARG_WRONG,… in ProcessOptions()
1340 …PetscCheck(!ctx->jacobian_field_major_order, ctx->comm, PETSC_ERR_ARG_WRONG, "-dm_landau_jacobian_… in ProcessOptions()
1343 …for (ii = ctx->num_species; ii < LANDAU_MAX_SPECIES; ii++) ctx->masses[ii] = ctx->thermal_temps[ii… in ProcessOptions()
1344 if (ctx->verbose != 0) { in ProcessOptions()
1347 … %10.3e %10.3e ...\n", (double)ctx->masses[0], (double)(ctx->masses[1] / pmassunit), (double)(ctx… in ProcessOptions()
1348 …10.3e\n", (double)ctx->charges[0], (double)(-ctx->charges[1] / ctx->charges[0]), (double)(ctx->num… in ProcessOptions()
1349 … i: %10.3e %10.3e\n", (double)ctx->n[0], (double)ctx->n[1], (double)(ctx->num_specie… in ProcessOptions()
1350 …0=%10.3e %" PetscInt_FMT " batched, view batch %" PetscInt_FMT "\n", (double)ctx->thermal_temps[0], in ProcessOptions()
1351 …ctx->thermal_temps[1], (double)((ctx->num_species > 2) ? ctx->thermal_temps[2] : 0), non_dim_grid,… in ProcessOptions()
1352 ….3e perp=%10.3e (%" PetscInt_FMT ") ", 0, (double)ctx->radius_par[0], (double)ctx->radius_perp[0],… in ProcessOptions()
1353 …ctx->num_grids; ii++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, ", %" PetscInt_FMT ": par=%10.3e per… in ProcessOptions()
1354 …if (ctx->use_relativistic_corrections) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nUse relativistic… in ProcessOptions()
1361 ctx->stage = 0; in ProcessOptions()
1362 PetscCall(PetscLogEventRegister("Landau Create", DM_CLASSID, &ctx->events[13])); /* 13 */ in ProcessOptions()
1363 PetscCall(PetscLogEventRegister(" GPU ass. setup", DM_CLASSID, &ctx->events[2])); /* 2 */ in ProcessOptions()
1364 PetscCall(PetscLogEventRegister(" Build matrix", DM_CLASSID, &ctx->events[12])); /* 12 */ in ProcessOptions()
1365 PetscCall(PetscLogEventRegister(" Assembly maps", DM_CLASSID, &ctx->events[15])); /* 15 */ in ProcessOptions()
1366 PetscCall(PetscLogEventRegister("Landau Mass mat", DM_CLASSID, &ctx->events[14])); /* 14 */ in ProcessOptions()
1367 PetscCall(PetscLogEventRegister("Landau Operator", DM_CLASSID, &ctx->events[11])); /* 11 */ in ProcessOptions()
1368 PetscCall(PetscLogEventRegister("Landau Jacobian", DM_CLASSID, &ctx->events[0])); /* 0 */ in ProcessOptions()
1369 PetscCall(PetscLogEventRegister("Landau Mass", DM_CLASSID, &ctx->events[9])); /* 9 */ in ProcessOptions()
1370 PetscCall(PetscLogEventRegister(" Preamble", DM_CLASSID, &ctx->events[10])); /* 10 */ in ProcessOptions()
1371 PetscCall(PetscLogEventRegister(" static IP Data", DM_CLASSID, &ctx->events[7])); /* 7 */ in ProcessOptions()
1372 PetscCall(PetscLogEventRegister(" dynamic IP-Jac", DM_CLASSID, &ctx->events[1])); /* 1 */ in ProcessOptions()
1373 PetscCall(PetscLogEventRegister(" Kernel-init", DM_CLASSID, &ctx->events[3])); /* 3 */ in ProcessOptions()
1374 PetscCall(PetscLogEventRegister(" Jac-f-df (GPU)", DM_CLASSID, &ctx->events[8])); /* 8 */ in ProcessOptions()
1375 PetscCall(PetscLogEventRegister(" J Kernel (GPU)", DM_CLASSID, &ctx->events[4])); /* 4 */ in ProcessOptions()
1376 PetscCall(PetscLogEventRegister(" M Kernel (GPU)", DM_CLASSID, &ctx->events[16])); /* 16 */ in ProcessOptions()
1377 PetscCall(PetscLogEventRegister(" Copy to CPU", DM_CLASSID, &ctx->events[5])); /* 5 */ in ProcessOptions()
1378 PetscCall(PetscLogEventRegister(" CPU assemble", DM_CLASSID, &ctx->events[6])); /* 6 */ in ProcessOptions()
1403 …rorCode CreateStaticData(PetscInt dim, IS grid_batch_is_inv[], const char prefix[], LandauCtx *ctx) in CreateStaticData() argument
1414 PetscCall(PetscFEGetDimension(ctx->fe[0], &Nb)); in CreateStaticData()
1415 …PetscCheck(Nb <= LANDAU_MAX_NQND, ctx->comm, PETSC_ERR_ARG_WRONG, "Order too high. Nb = %" PetscIn… in CreateStaticData()
1416 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in CreateStaticData()
1417 for (PetscInt ii = ctx->species_offset[grid]; ii < ctx->species_offset[grid + 1]; ii++) { in CreateStaticData()
1418 invMass[ii] = ctx->m_0 / ctx->masses[ii]; in CreateStaticData()
1419 nu_alpha[ii] = PetscSqr(ctx->charges[ii] / ctx->m_0) * ctx->m_0 / ctx->masses[ii]; in CreateStaticData()
1420 …nu_beta[ii] = PetscSqr(ctx->charges[ii] / ctx->epsilon0) / (8 * PETSC_PI) * ctx->t_0 * ctx->n_0 /… in CreateStaticData()
1423 if (ctx->verbose == 4) { in CreateStaticData()
1425 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in CreateStaticData()
1426 PetscInt iii = ctx->species_offset[grid]; in CreateStaticData()
1427 …for (PetscInt ii = iii; ii < ctx->species_offset[grid + 1]; ii++) PetscCall(PetscPrintf(PETSC_COMM… in CreateStaticData()
1430 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in CreateStaticData()
1431 PetscInt iii = ctx->species_offset[grid]; in CreateStaticData()
1432 …for (PetscInt ii = iii; ii < ctx->species_offset[grid + 1]; ii++) PetscCall(PetscPrintf(PETSC_COMM… in CreateStaticData()
1435 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in CreateStaticData()
1436 PetscInt iii = ctx->species_offset[grid]; in CreateStaticData()
1437 for (PetscInt ii = iii; ii < ctx->species_offset[grid + 1]; ii++) { in CreateStaticData()
1438 for (PetscInt gridj = 0; gridj < ctx->num_grids; gridj++) { in CreateStaticData()
1439 PetscInt jjj = ctx->species_offset[gridj]; in CreateStaticData()
1440 … jjj; jj < ctx->species_offset[gridj + 1]; jj++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %14.9e"… in CreateStaticData()
1446 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in CreateStaticData()
1447 PetscInt iii = ctx->species_offset[grid]; in CreateStaticData()
1448 for (PetscInt ii = iii; ii < ctx->species_offset[grid + 1]; ii++) { in CreateStaticData()
1449 for (PetscInt gridj = 0; gridj < ctx->num_grids; gridj++) { in CreateStaticData()
1450 PetscInt jjj = ctx->species_offset[gridj]; in CreateStaticData()
1451 … (PetscInt jj = jjj; jj < ctx->species_offset[gridj + 1]; jj++) PetscCall(PetscPrintf(PETSC_COMM_W… in CreateStaticData()
1457 PetscCall(DMGetDS(ctx->plex[0], &prob)); // same DS for all grids in CreateStaticData()
1460 PetscCheck(ctx->plex[0], ctx->comm, PETSC_ERR_ARG_WRONG, "Plex not created"); in CreateStaticData()
1461 PetscCall(PetscFEGetQuadrature(ctx->fe[0], &quad)); in CreateStaticData()
1463 …PetscCheck(Nq <= LANDAU_MAX_NQND, ctx->comm, PETSC_ERR_ARG_WRONG, "Order too high. Nq = %" PetscIn… in CreateStaticData()
1465 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in CreateStaticData()
1467 PetscCheck(ctx->plex[grid] != NULL, ctx->comm, PETSC_ERR_ARG_WRONG, "Plex not created"); in CreateStaticData()
1468 PetscCall(DMPlexGetHeightStratum(ctx->plex[grid], 0, &cStart, &cEnd)); in CreateStaticData()
1470 PetscCall(DMGetLocalSection(ctx->plex[grid], §ion[grid])); in CreateStaticData()
1471 PetscCall(DMGetGlobalSection(ctx->plex[grid], &globsection[grid])); in CreateStaticData()
1476 if (ctx->gpu_assembly) { /* we need GPU object with GPU assembly */ in CreateStaticData()
1481 const PetscInt *plex_batch = NULL, elMatSz = Nb * Nb * ctx->num_species * ctx->num_species; in CreateStaticData()
1484 PetscCall(PetscInfo(ctx->plex[0], "Make GPU maps %d\n", 1)); in CreateStaticData()
1485 PetscCall(PetscLogEventBegin(ctx->events[2], 0, 0, 0, 0)); in CreateStaticData()
1486 PetscCall(PetscMalloc(sizeof(*maps) * ctx->num_grids, &maps)); in CreateStaticData()
1493 … PetscCall(PetscInfo(ctx->plex[0], "COO initialization, %" PetscInt_FMT " cells\n", ncellsTot)); in CreateStaticData()
1494 ctx->SData_d.coo_n_cellsTot = ncellsTot; in CreateStaticData()
1495 ctx->SData_d.coo_elem_offsets = (void *)coo_elem_offsets; in CreateStaticData()
1496 ctx->SData_d.coo_elem_fullNb = (void *)coo_elem_fullNb; in CreateStaticData()
1497 ctx->SData_d.coo_elem_point_offsets = (void *)coo_elem_point_offsets; in CreateStaticData()
1500 ctx->SData_d.coo_max_fullnb = 0; in CreateStaticData()
1501 for (PetscInt grid = 0, glb_elem_idx = 0; grid < ctx->num_grids; grid++) { in CreateStaticData()
1504 …PetscCheck(!plex_batch, ctx->comm, PETSC_ERR_ARG_WRONG, "-dm_landau_jacobian_field_major_order DEP… in CreateStaticData()
1505 PetscCall(DMPlexGetHeightStratum(ctx->plex[grid], 0, &cStart, &cEnd)); in CreateStaticData()
1512 maps[grid].deviceType = ctx->deviceType; in CreateStaticData()
1513 maps[grid].numgrids = ctx->num_grids; in CreateStaticData()
1525 …PetscCall(DMPlexGetClosureIndices(ctx->plex[grid], section[grid], globsection[grid], ej, PETSC_TRU… in CreateStaticData()
1526 if (ctx->simplex) { in CreateStaticData()
1527 …PetscCheck(numindices == Nb, ctx->comm, PETSC_ERR_ARG_WRONG, "numindices != Nb numindices=%" Petsc… in CreateStaticData()
1546 … PetscCheck(!ctx->simplex, ctx->comm, PETSC_ERR_ARG_WRONG, "No constraints with simplex"); in CreateStaticData()
1579 … if (tmp != 0) PetscCall(PetscPrintf(ctx->comm, " | %22.16e\n", (double)tmp)); in CreateStaticData()
1590 …PetscCall(DMPlexRestoreClosureIndices(ctx->plex[grid], section[grid], globsection[grid], ej, PETSC… in CreateStaticData()
1591 …if (elMat != valuesOrig) PetscCall(DMRestoreWorkArray(ctx->plex[grid], numindices * numindices, MP… in CreateStaticData()
1597 if (fullNb > ctx->SData_d.coo_max_fullnb) ctx->SData_d.coo_max_fullnb = fullNb; in CreateStaticData()
1611 …if (ctx->deviceType == LANDAU_KOKKOS) PetscCall(LandauKokkosCreateMatMaps(maps, pointMaps, Nf, gri… in CreateStaticData()
1621 ctx->SData_d.coo_size = coo_elem_offsets[ncellsTot] * ctx->batch_sz; in CreateStaticData()
1622 PetscCall(PetscMalloc2(ctx->SData_d.coo_size, &oor, ctx->SData_d.coo_size, &ooc)); in CreateStaticData()
1623 for (PetscInt i = 0; i < ctx->SData_d.coo_size; i++) oor[i] = ooc[i] = -1; in CreateStaticData()
1625 for (PetscInt grid = 0, glb_elem_idx = 0; grid < ctx->num_grids; grid++) { in CreateStaticData()
1650 for (PetscInt b_id = 0; b_id < ctx->batch_sz; b_id++) { in CreateStaticData()
1651 for (PetscInt grid = 0, glb_elem_idx = 0; grid < ctx->num_grids; grid++) { in CreateStaticData()
1652 … const PetscInt moffset = LAND_MOFFSET(b_id, grid, ctx->batch_sz, ctx->num_grids, ctx->mat_offset); in CreateStaticData()
1686 PetscCall(MatSetPreallocationCOO(ctx->J, ctx->SData_d.coo_size, oor, ooc)); in CreateStaticData()
1694 PetscCall(PetscObjectCompose((PetscObject)ctx->J, "assembly_maps", (PetscObject)container)); in CreateStaticData()
1696 PetscCall(PetscLogEventEnd(ctx->events[2], 0, 0, 0, 0)); in CreateStaticData()
1702 PetscCall(PetscLogEventBegin(ctx->events[7], 0, 0, 0, 0)); in CreateStaticData()
1703 PetscCall(PetscInfo(ctx->plex[0], "Initialize static data\n")); in CreateStaticData()
1704 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) nip_glb += Nq * numCells[grid]; in CreateStaticData()
1706 if (ctx->verbose != 0) { in CreateStaticData()
1709 PetscCall(MatGetInfo(ctx->J, MAT_LOCAL, &info)); in CreateStaticData()
1710 PetscCall(MatGetSize(ctx->J, &N, NULL)); in CreateStaticData()
1711 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) ncells += numCells[grid]; in CreateStaticData()
1713 ctx->num_species, Nb, dim, N, (PetscInt)info.nz_used)); in CreateStaticData()
1717 if (ctx->use_energy_tensor_trick) { in CreateStaticData()
1718 … PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, ctx->simplex, prefix, PETSC_DECIDE, &fe)); in CreateStaticData()
1722 for (grid = 0, outer_ipidx = 0, outer_ej = 0; grid < ctx->num_grids; grid++) { // OpenMP (once) in CreateStaticData()
1728 PetscCall(DMPlexGetHeightStratum(ctx->plex[grid], 0, &cStart, &cEnd)); in CreateStaticData()
1730 if (ctx->use_energy_tensor_trick) { in CreateStaticData()
1731 …etscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *) = {ctx->use_relativistic_c… in CreateStaticData()
1733 PetscReal *c2_0[1], data[1] = {PetscSqr(C_0(ctx->v_0))}; in CreateStaticData()
1735 PetscCall(DMClone(ctx->plex[grid], &dmEnergy)); in CreateStaticData()
1755 …DAU_MAX_NQND], Jdummy[LANDAU_MAX_NQND * LANDAU_DIM * LANDAU_DIM], c0 = C_0(ctx->v_0), c02 = PetscS… in CreateStaticData()
1757 …PetscCall(DMPlexComputeCellGeometryFEM(ctx->plex[grid], ej + cStart, quad, vj, Jdummy, invJe, detJ… in CreateStaticData()
1758 …if (ctx->use_energy_tensor_trick) PetscCall(DMPlexVecGetClosure(dmEnergy, e_section, v2_2, ej + cS… in CreateStaticData()
1766 if (ctx->use_energy_tensor_trick) { in CreateStaticData()
1775 if (ctx->use_relativistic_corrections) { in CreateStaticData()
1784 …PetscCall(PetscPrintf(ctx->comm, "Error: %12.5e %" PetscInt_FMT ".%" PetscInt_FMT ") dg2/c02 = %12… in CreateStaticData()
1805 …if (ctx->use_energy_tensor_trick) PetscCall(DMPlexVecRestoreClosure(dmEnergy, e_section, v2_2, ej … in CreateStaticData()
1807 if (ctx->use_energy_tensor_trick) { in CreateStaticData()
1812 if (ctx->use_energy_tensor_trick) PetscCall(PetscFEDestroy(&fe)); in CreateStaticData()
1814 if (ctx->deviceType == LANDAU_KOKKOS) { in CreateStaticData()
1816 …ctx->plex[0], Nq, Nb, ctx->batch_sz, ctx->num_grids, numCells, ctx->species_offset, ctx->mat_offse… in CreateStaticData()
1821 SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONG, "-landau_device_type kokkos not built"); in CreateStaticData()
1824 …l *nu_alpha_p = (PetscReal *)ctx->SData_d.alpha, *nu_beta_p = (PetscReal *)ctx->SData_d.beta, *inv… in CreateStaticData()
1825 ctx->SData_d.w = (void *)ww; in CreateStaticData()
1826 ctx->SData_d.x = (void *)xx; in CreateStaticData()
1827 ctx->SData_d.y = (void *)yy; in CreateStaticData()
1828 ctx->SData_d.z = (void *)zz; in CreateStaticData()
1829 ctx->SData_d.invJ = (void *)invJ_a; in CreateStaticData()
1830 …PetscCall(PetscMalloc4(ctx->num_species, &nu_alpha_p, ctx->num_species, &nu_beta_p, ctx->num_speci… in CreateStaticData()
1831 for (PetscInt ii = 0; ii < ctx->num_species; ii++) { in CreateStaticData()
1836 ctx->SData_d.alpha = (void *)nu_alpha_p; in CreateStaticData()
1837 ctx->SData_d.beta = (void *)nu_beta_p; in CreateStaticData()
1838 ctx->SData_d.invMass = (void *)invMass_p; in CreateStaticData()
1839 ctx->SData_d.lambdas = (void *)lambdas_p; in CreateStaticData()
1841 …GRIDS][LANDAU_MAX_GRIDS] = (PetscReal (*)[LANDAU_MAX_GRIDS][LANDAU_MAX_GRIDS])ctx->SData_d.lambdas; in CreateStaticData()
1842 …for (PetscInt gridj = 0; gridj < LANDAU_MAX_GRIDS; gridj++) (*lambdas)[grid][gridj] = ctx->lambdas… in CreateStaticData()
1845 PetscCall(PetscLogEventEnd(ctx->events[7], 0, 0, 0, 0)); in CreateStaticData()
1874 …uCreateJacobianMatrix(MPI_Comm comm, Vec X, IS grid_batch_is_inv[LANDAU_MAX_GRIDS], LandauCtx *ctx) in LandauCreateJacobianMatrix() argument
1880 if (!ctx->gpu_assembly) { /* we need GPU object with GPU assembly */ in LandauCreateJacobianMatrix()
1884 …if (ctx->gpu_assembly && ctx->jacobian_field_major_order) PetscCall(PetscMalloc1(ctx->mat_offset[c… in LandauCreateJacobianMatrix()
1885 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauCreateJacobianMatrix()
1886 const PetscInt *values, n = ctx->mat_offset[grid + 1] - ctx->mat_offset[grid]; in LandauCreateJacobianMatrix()
1892 PetscCall(DMClone(ctx->plex[grid], &massDM)); in LandauCreateJacobianMatrix()
1893 PetscCall(DMCopyFields(ctx->plex[grid], PETSC_DETERMINE, PETSC_DETERMINE, massDM)); in LandauCreateJacobianMatrix()
1896 …for (PetscInt ix = 0, ii = ctx->species_offset[grid]; ii < ctx->species_offset[grid + 1]; ii++, ix… in LandauCreateJacobianMatrix()
1902 PetscCall(DMCreateLocalVector(ctx->plex[grid], &tvec)); in LandauCreateJacobianMatrix()
1903 PetscCall(DMPlexSNESComputeJacobianFEM(massDM, tvec, gMat, gMat, ctx)); in LandauCreateJacobianMatrix()
1908 if (ctx->gpu_assembly && ctx->jacobian_field_major_order) { in LandauCreateJacobianMatrix()
1914 …for (PetscInt b_id = 0; b_id < ctx->batch_sz; b_id++) { // add batch size DMs for this species grid in LandauCreateJacobianMatrix()
1916 PetscInt N = ctx->mat_offset[ctx->num_grids], n0 = ctx->mat_offset[grid] + b_id * N; in LandauCreateJacobianMatrix()
1919 PetscInt n0 = ctx->mat_offset[grid] * ctx->batch_sz + b_id * n; in LandauCreateJacobianMatrix()
1928 …ctx->gpu_assembly && ctx->jacobian_field_major_order) PetscCall(ISCreateGeneral(comm, ctx->mat_off… in LandauCreateJacobianMatrix()
1930 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauCreateJacobianMatrix()
1935 for (PetscInt b_id = 0; b_id < ctx->batch_sz; b_id++) { in LandauCreateJacobianMatrix()
1936 …const PetscInt moffset = LAND_MOFFSET(b_id, grid, ctx->batch_sz, ctx->num_grids, ctx->mat_offs… in LandauCreateJacobianMatrix()
1943 …PetscCall(PetscInfo(ctx->plex[grid], "Realloc buffer %" PetscInt_FMT " to %" PetscInt_FMT " (row s… in LandauCreateJacobianMatrix()
1950 PetscCall(MatSetValues(ctx->J, 1, &row, nzl, colbuf, vals, INSERT_VALUES)); in LandauCreateJacobianMatrix()
1956 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) PetscCall(MatDestroy(&subM[grid])); in LandauCreateJacobianMatrix()
1957 PetscCall(MatAssemblyBegin(ctx->J, MAT_FINAL_ASSEMBLY)); in LandauCreateJacobianMatrix()
1958 PetscCall(MatAssemblyEnd(ctx->J, MAT_FINAL_ASSEMBLY)); in LandauCreateJacobianMatrix()
1961 PetscCall(MatViewFromOptions(ctx->J, NULL, "-dm_landau_mat_view")); in LandauCreateJacobianMatrix()
1962 if (ctx->gpu_assembly && ctx->jacobian_field_major_order) { in LandauCreateJacobianMatrix()
1964 …PetscCall(MatCreateSubMatrix(ctx->J, ctx->batch_is, ctx->batch_is, MAT_INITIAL_MATRIX, &mat_block_… in LandauCreateJacobianMatrix()
1967 PetscCall(VecScatterCreate(X, ctx->batch_is, X, NULL, &ctx->plex_batch)); in LandauCreateJacobianMatrix()
1968 PetscCall(VecDuplicate(X, &ctx->work_vec)); in LandauCreateJacobianMatrix()
2051 LandauCtx *ctx; in DMPlexLandauCreateVelocitySpace() local
2058 PetscCall(PetscNew(&ctx)); in DMPlexLandauCreateVelocitySpace()
2059 ctx->comm = comm; /* used for diagnostics and global errors */ in DMPlexLandauCreateVelocitySpace()
2061 PetscCall(ProcessOptions(ctx, prefix)); in DMPlexLandauCreateVelocitySpace()
2062 if (dim == 2) ctx->use_relativistic_corrections = PETSC_FALSE; in DMPlexLandauCreateVelocitySpace()
2065 PetscCall(PetscLogEventBegin(ctx->events[13], 0, 0, 0, 0)); in DMPlexLandauCreateVelocitySpace()
2066 PetscCall(PetscLogEventBegin(ctx->events[15], 0, 0, 0, 0)); in DMPlexLandauCreateVelocitySpace()
2067 …PetscCall(LandauDMCreateVMeshes(PETSC_COMM_SELF, dim, prefix, ctx, *pack)); // creates grids (Fore… in DMPlexLandauCreateVelocitySpace()
2068 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauCreateVelocitySpace()
2070 PetscCall(SetupDS(ctx->plex[grid], dim, grid, prefix, ctx)); in DMPlexLandauCreateVelocitySpace()
2072 PetscCall(DMCreateGlobalVector(ctx->plex[grid], &Xsub[grid])); in DMPlexLandauCreateVelocitySpace()
2075 PetscCall(LandauSetInitialCondition(ctx->plex[grid], Xsub[grid], grid, 0, 1, ctx)); in DMPlexLandauCreateVelocitySpace()
2077 if (ctx->use_p4est) { in DMPlexLandauCreateVelocitySpace()
2079 PetscCall(adapt(grid, ctx, &Xsub[grid])); // forest goes in, plex comes out in DMPlexLandauCreateVelocitySpace()
2081 PetscCall(DMConvert(ctx->plex[grid], DMPLEX, &plex)); in DMPlexLandauCreateVelocitySpace()
2082 PetscCall(DMDestroy(&ctx->plex[grid])); in DMPlexLandauCreateVelocitySpace()
2083 ctx->plex[grid] = plex; in DMPlexLandauCreateVelocitySpace()
2084 } else if (ctx->sphere && dim == 3) { in DMPlexLandauCreateVelocitySpace()
2085 …if (ctx->map_sphere) PetscCall(LandauSphereMesh(ctx->plex[grid], ctx->radius[grid] * ctx->sphere_i… in DMPlexLandauCreateVelocitySpace()
2086 PetscCall(LandauSetInitialCondition(ctx->plex[grid], Xsub[grid], grid, 0, 1, ctx)); in DMPlexLandauCreateVelocitySpace()
2089 PetscCall(DMViewFromOptions(ctx->plex[grid], NULL, "-dm_landau_amr_dm_view")); in DMPlexLandauCreateVelocitySpace()
2094 PetscCall(DMCompositeAddDM(*pack, ctx->plex[grid])); in DMPlexLandauCreateVelocitySpace()
2096 …for (PetscInt b_id = 0; b_id < ctx->batch_sz; b_id++) { // add batch size DMs for this species grid in DMPlexLandauCreateVelocitySpace()
2097 PetscCall(DMCompositeAddDM(*pack, ctx->plex[grid])); in DMPlexLandauCreateVelocitySpace()
2100 PetscCall(DMSetApplicationContext(ctx->plex[grid], ctx)); in DMPlexLandauCreateVelocitySpace()
2104 for (PetscInt b_id = 1; b_id < ctx->batch_sz; b_id++) { in DMPlexLandauCreateVelocitySpace()
2105 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) PetscCall(DMCompositeAddDM(*pack, ctx->plex… in DMPlexLandauCreateVelocitySpace()
2109 ctx->mat_offset[0] = 0; in DMPlexLandauCreateVelocitySpace()
2110 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauCreateVelocitySpace()
2113 ctx->mat_offset[grid + 1] = ctx->mat_offset[grid] + n; in DMPlexLandauCreateVelocitySpace()
2116 PetscCall(DMSetApplicationContext(*pack, ctx)); in DMPlexLandauCreateVelocitySpace()
2118 PetscCall(DMCreateMatrix(*pack, &ctx->J)); in DMPlexLandauCreateVelocitySpace()
2120 PetscCall(MatSetOption(ctx->J, MAT_STRUCTURALLY_SYMMETRIC, PETSC_TRUE)); in DMPlexLandauCreateVelocitySpace()
2121 PetscCall(MatSetOption(ctx->J, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE)); in DMPlexLandauCreateVelocitySpace()
2122 PetscCall(PetscObjectSetName((PetscObject)ctx->J, "Jac")); in DMPlexLandauCreateVelocitySpace()
2125 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauCreateVelocitySpace()
2128 for (PetscInt b_id = 0; b_id < ctx->batch_sz; b_id++) { in DMPlexLandauCreateVelocitySpace()
2130 …const PetscInt moffset = LAND_MOFFSET(b_id, grid, ctx->batch_sz, ctx->num_grids, ctx->mat_offs… in DMPlexLandauCreateVelocitySpace()
2131 … PetscCall(LandauSetInitialCondition(ctx->plex[grid], Xsub[grid], grid, b_id, ctx->batch_sz, ctx)); in DMPlexLandauCreateVelocitySpace()
2138 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) PetscCall(VecDestroy(&Xsub[grid])); in DMPlexLandauCreateVelocitySpace()
2140 if (ctx->gpu_assembly) { /* we need GPU object with GPU assembly */ in DMPlexLandauCreateVelocitySpace()
2142 if (ctx->deviceType == LANDAU_KOKKOS) { in DMPlexLandauCreateVelocitySpace()
2143 …PetscCall(PetscObjectTypeCompareAny((PetscObject)ctx->J, &flg, MATSEQAIJKOKKOS, MATMPIAIJKOKKOS, M… in DMPlexLandauCreateVelocitySpace()
2145 …PetscCheck(flg, ctx->comm, PETSC_ERR_ARG_WRONG, "must use '-dm_mat_type aijkokkos -dm_vec_type kok… in DMPlexLandauCreateVelocitySpace()
2147 …PetscCheck(flg, ctx->comm, PETSC_ERR_ARG_WRONG, "must configure with '--download-kokkos-kernels' f… in DMPlexLandauCreateVelocitySpace()
2151 PetscCall(PetscLogEventEnd(ctx->events[15], 0, 0, 0, 0)); in DMPlexLandauCreateVelocitySpace()
2154 ctx->work_vec = NULL; in DMPlexLandauCreateVelocitySpace()
2155 ctx->plex_batch = NULL; in DMPlexLandauCreateVelocitySpace()
2156 ctx->batch_is = NULL; in DMPlexLandauCreateVelocitySpace()
2158 PetscCall(PetscLogEventBegin(ctx->events[12], 0, 0, 0, 0)); in DMPlexLandauCreateVelocitySpace()
2159 PetscCall(LandauCreateJacobianMatrix(comm, *X, grid_batch_is_inv, ctx)); in DMPlexLandauCreateVelocitySpace()
2160 PetscCall(PetscLogEventEnd(ctx->events[12], 0, 0, 0, 0)); in DMPlexLandauCreateVelocitySpace()
2163 PetscCall(CreateStaticData(dim, grid_batch_is_inv, prefix, ctx)); in DMPlexLandauCreateVelocitySpace()
2165 PetscCall(PetscLogEventEnd(ctx->events[13], 0, 0, 0, 0)); in DMPlexLandauCreateVelocitySpace()
2170 if (J) *J = ctx->J; in DMPlexLandauCreateVelocitySpace()
2172 if (ctx->gpu_assembly && ctx->jacobian_field_major_order) { in DMPlexLandauCreateVelocitySpace()
2176 PetscCall(PetscContainerSetPointer(container, (void *)ctx)); in DMPlexLandauCreateVelocitySpace()
2177 PetscCall(PetscObjectCompose((PetscObject)ctx->J, "LandauCtx", (PetscObject)container)); in DMPlexLandauCreateVelocitySpace()
2181 PetscCall(PetscContainerSetPointer(container, (void *)ctx->plex_batch)); in DMPlexLandauCreateVelocitySpace()
2182 PetscCall(PetscObjectCompose((PetscObject)ctx->J, "plex_batch_is", (PetscObject)container)); in DMPlexLandauCreateVelocitySpace()
2191 *pNf = ctx->batch_sz; in DMPlexLandauCreateVelocitySpace()
2194 PetscCall(PetscObjectCompose((PetscObject)ctx->J, "batch size", (PetscObject)container)); in DMPlexLandauCreateVelocitySpace()
2219 LandauCtx *ctx; in DMPlexLandauAccess() local
2222 …PetscCall(DMGetApplicationContext(pack, &ctx)); // uses ctx->num_grids; ctx->plex[grid]; ctx->batc… in DMPlexLandauAccess()
2223 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauAccess()
2226 …for (PetscInt sp = ctx->species_offset[grid], i0 = 0; sp < ctx->species_offset[grid + 1]; sp++, i0… in DMPlexLandauAccess()
2231 PetscCall(DMCreateSubDM(ctx->plex[grid], 1, vf, &vis, &vdm)); in DMPlexLandauAccess()
2232 PetscCall(DMSetApplicationContext(vdm, ctx)); // the user might want this in DMPlexLandauAccess()
2235 for (PetscInt b_id = 0; b_id < ctx->batch_sz; b_id++) { in DMPlexLandauAccess()
2236 … const PetscInt moffset = LAND_MOFFSET(b_id, grid, ctx->batch_sz, ctx->num_grids, ctx->mat_offset); in DMPlexLandauAccess()
2271 LandauCtx *ctx; in DMPlexLandauDestroyVelocitySpace() local
2274 PetscCall(DMGetApplicationContext(*dm, &ctx)); in DMPlexLandauDestroyVelocitySpace()
2275 PetscCall(MatDestroy(&ctx->M)); in DMPlexLandauDestroyVelocitySpace()
2276 PetscCall(MatDestroy(&ctx->J)); in DMPlexLandauDestroyVelocitySpace()
2277 for (PetscInt ii = 0; ii < ctx->num_species; ii++) PetscCall(PetscFEDestroy(&ctx->fe[ii])); in DMPlexLandauDestroyVelocitySpace()
2278 PetscCall(ISDestroy(&ctx->batch_is)); in DMPlexLandauDestroyVelocitySpace()
2279 PetscCall(VecDestroy(&ctx->work_vec)); in DMPlexLandauDestroyVelocitySpace()
2280 PetscCall(VecScatterDestroy(&ctx->plex_batch)); in DMPlexLandauDestroyVelocitySpace()
2281 if (ctx->deviceType == LANDAU_KOKKOS) { in DMPlexLandauDestroyVelocitySpace()
2283 PetscCall(LandauKokkosStaticDataClear(&ctx->SData_d)); in DMPlexLandauDestroyVelocitySpace()
2285 SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONG, "-landau_device_type %s not built", "kokkos"); in DMPlexLandauDestroyVelocitySpace()
2288 if (ctx->SData_d.x) { /* in a CPU run */ in DMPlexLandauDestroyVelocitySpace()
2289 … *)ctx->SData_d.invJ, *xx = (PetscReal *)ctx->SData_d.x, *yy = (PetscReal *)ctx->SData_d.y, *zz = … in DMPlexLandauDestroyVelocitySpace()
2290 …ctx->SData_d.coo_elem_offsets, *coo_elem_fullNb = (LandauIdx *)ctx->SData_d.coo_elem_fullNb, (*coo… in DMPlexLandauDestroyVelocitySpace()
2294 …PetscCall(PetscFree4(ctx->SData_d.alpha, ctx->SData_d.beta, ctx->SData_d.invMass, ctx->SData_d.lam… in DMPlexLandauDestroyVelocitySpace()
2298 if (ctx->times[LANDAU_MATRIX_TOTAL] > 0) { // OMP timings in DMPlexLandauDestroyVelocitySpace()
2299 …PetscCall(PetscPrintf(ctx->comm, "TSStep N 1.0 %10.3e\n", ctx->times[LANDAU_EX2_TSS… in DMPlexLandauDestroyVelocitySpace()
2300 …ntf(ctx->comm, "2: Solve: %10.3e with %" PetscInt_FMT " threads\n", ctx->times[LANDAU_E… in DMPlexLandauDestroyVelocitySpace()
2301 …PetscCall(PetscPrintf(ctx->comm, "3: Landau: %10.3e\n", ctx->times[LANDAU_MATRIX_TOTAL])… in DMPlexLandauDestroyVelocitySpace()
2302 …scCall(PetscPrintf(ctx->comm, "Landau Jacobian %" PetscInt_FMT " 1.0 %10.3e\n", (PetscInt)ct… in DMPlexLandauDestroyVelocitySpace()
2303 …PetscCall(PetscPrintf(ctx->comm, "Landau Operator N 1.0 %10.3e\n", ctx->times[LANDAU_OPERAT… in DMPlexLandauDestroyVelocitySpace()
2304 …PetscCall(PetscPrintf(ctx->comm, "Landau Mass N 1.0 %10.3e\n", ctx->times[LANDAU_MASS])… in DMPlexLandauDestroyVelocitySpace()
2305 …PetscCall(PetscPrintf(ctx->comm, " Jac-f-df (GPU) N 1.0 %10.3e\n", ctx->times[LANDAU_F_DF])… in DMPlexLandauDestroyVelocitySpace()
2306 …PetscCall(PetscPrintf(ctx->comm, " Kernel (GPU) N 1.0 %10.3e\n", ctx->times[LANDAU_KERNEL… in DMPlexLandauDestroyVelocitySpace()
2307 … PetscCall(PetscPrintf(ctx->comm, "MatLUFactorNum X 1.0 %10.3e\n", ctx->times[KSP_FACTOR])); in DMPlexLandauDestroyVelocitySpace()
2308 … PetscCall(PetscPrintf(ctx->comm, "MatSolve X 1.0 %10.3e\n", ctx->times[KSP_SOLVE])); in DMPlexLandauDestroyVelocitySpace()
2310 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) PetscCall(DMDestroy(&ctx->plex[grid])); in DMPlexLandauDestroyVelocitySpace()
2311 PetscCall(PetscFree(ctx)); in DMPlexLandauDestroyVelocitySpace()
2394 LandauCtx *ctx; in DMPlexLandauPrintNorms() local
2407 PetscCall(DMGetApplicationContext(pack, &ctx)); in DMPlexLandauPrintNorms()
2408 PetscCheck(ctx, PETSC_COMM_SELF, PETSC_ERR_PLIB, "no context"); in DMPlexLandauPrintNorms()
2411 …DMs == ctx->num_grids * ctx->batch_sz, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "#DM wrong %" PetscInt_FM… in DMPlexLandauPrintNorms()
2414 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauPrintNorms()
2415 Vec Xloc = globXArray[LAND_PACK_IDX(ctx->batch_view_idx, grid)]; in DMPlexLandauPrintNorms()
2416 PetscCall(DMGetDS(ctx->plex[grid], &prob)); in DMPlexLandauPrintNorms()
2417 for (ii = ctx->species_offset[grid], i0 = 0; ii < ctx->species_offset[grid + 1]; ii++, i0++) { in DMPlexLandauPrintNorms()
2418 PetscScalar user[2] = {(PetscScalar)i0, ctx->charges[ii]}; in DMPlexLandauPrintNorms()
2422 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2423 density[ii] = tt[0] * ctx->n_0 * ctx->charges[ii]; in DMPlexLandauPrintNorms()
2425 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2426 zmomentum[ii] = tt[0] * ctx->n_0 * ctx->v_0 * ctx->masses[ii]; in DMPlexLandauPrintNorms()
2428 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2429 energy[ii] = tt[0] * 0.5 * ctx->n_0 * ctx->v_0 * ctx->v_0 * ctx->masses[ii]; in DMPlexLandauPrintNorms()
2436 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2437 density[ii] = tt[0] * ctx->n_0 * ctx->charges[ii]; in DMPlexLandauPrintNorms()
2441 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2442 xmomentum[ii] = tt[0] * ctx->n_0 * ctx->v_0 * ctx->masses[ii]; in DMPlexLandauPrintNorms()
2445 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2446 ymomentum[ii] = tt[0] * ctx->n_0 * ctx->v_0 * ctx->masses[ii]; in DMPlexLandauPrintNorms()
2449 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2450 zmomentum[ii] = tt[0] * ctx->n_0 * ctx->v_0 * ctx->masses[ii]; in DMPlexLandauPrintNorms()
2451 if (ctx->use_relativistic_corrections) { in DMPlexLandauPrintNorms()
2463 PetscCall(MatMult(ctx->M, X, Mf)); in DMPlexLandauPrintNorms()
2466 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // yes a grid loop in a grid loop to prin… in DMPlexLandauPrintNorms()
2467 Vec v1 = globGammaArray[LAND_PACK_IDX(ctx->batch_view_idx, grid)]; in DMPlexLandauPrintNorms()
2468 data[0] = PetscSqr(C_0(ctx->v_0)); in DMPlexLandauPrintNorms()
2470 … PetscCall(DMProjectFunction(ctx->plex[grid], 0., gammaf, (void **)c2_0, INSERT_ALL_VALUES, v1)); in DMPlexLandauPrintNorms()
2476 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // yes a grid loop in a grid loop to prin… in DMPlexLandauPrintNorms()
2477 PetscInt Nf = ctx->species_offset[grid + 1] - ctx->species_offset[grid], N, bs; in DMPlexLandauPrintNorms()
2478 …ec Mfsub = globMfArray[LAND_PACK_IDX(ctx->batch_view_idx, grid)], Gsub = globGammaArray[LAND_… in DMPlexLandauPrintNorms()
2481 PetscCall(VecCreate(ctx->comm, &v1)); in DMPlexLandauPrintNorms()
2483 PetscCall(VecCreate(ctx->comm, &v2)); in DMPlexLandauPrintNorms()
2492 for (PetscInt i = 0, ix = ctx->species_offset[grid]; i < Nf; i++, ix++) { in DMPlexLandauPrintNorms()
2497 energy[ix] = PetscRealPart(val) * ctx->n_0 * ctx->v_0 * ctx->v_0 * ctx->masses[ix]; in DMPlexLandauPrintNorms()
2511 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2512 energy[ii] = 0.5 * tt[0] * ctx->n_0 * ctx->v_0 * ctx->v_0 * ctx->masses[ii]; in DMPlexLandauPrintNorms()
2521 if (ctx->num_species > 1) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n")); in DMPlexLandauPrintNorms()
2527 PetscCall(DMPlexGetHeightStratum(ctx->plex[0], 0, &cStart, &cEnd)); in DMPlexLandauPrintNorms()
2528 if (ctx->num_species > 1) { in DMPlexLandauPrintNorms()
2531 (double)(ctx->masses[1] / ctx->masses[0]), cEnd - cStart)); in DMPlexLandauPrintNorms()
2534 (double)(ctx->masses[1] / ctx->masses[0]), cEnd - cStart)); in DMPlexLandauPrintNorms()
2562 LandauCtx *ctx; in DMPlexLandauCreateMassMatrix() local
2568 PetscCall(DMGetApplicationContext(pack, &ctx)); in DMPlexLandauCreateMassMatrix()
2569 PetscCheck(ctx, PETSC_COMM_SELF, PETSC_ERR_PLIB, "no context"); in DMPlexLandauCreateMassMatrix()
2570 PetscCall(PetscLogEventBegin(ctx->events[14], 0, 0, 0, 0)); in DMPlexLandauCreateMassMatrix()
2574 for (PetscInt grid = 0, ix = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauCreateMassMatrix()
2575 PetscCall(DMClone(ctx->plex[grid], &massDM[grid])); in DMPlexLandauCreateMassMatrix()
2576 PetscCall(DMCopyFields(ctx->plex[grid], PETSC_DETERMINE, PETSC_DETERMINE, massDM[grid])); in DMPlexLandauCreateMassMatrix()
2579 for (ix = 0, ii = ctx->species_offset[grid]; ii < ctx->species_offset[grid + 1]; ii++, ix++) { in DMPlexLandauCreateMassMatrix()
2586 …for (PetscInt b_id = 0; b_id < ctx->batch_sz; b_id++) { // add batch size DMs for this species grid in DMPlexLandauCreateMassMatrix()
2594 for (PetscInt b_id = 1; b_id < ctx->batch_sz; b_id++) { in DMPlexLandauCreateMassMatrix()
2595 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) PetscCall(DMCompositeAddDM(mass_pack, massD… in DMPlexLandauCreateMassMatrix()
2605 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauCreateMassMatrix()
2610 PetscCall(DMPlexSNESComputeJacobianFEM(plex, locX, subM[grid], subM[grid], ctx)); in DMPlexLandauCreateMassMatrix()
2614 PetscCall(MatGetSize(ctx->J, &N1, NULL)); in DMPlexLandauCreateMassMatrix()
2618 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauCreateMassMatrix()
2623 for (PetscInt b_id = 0; b_id < ctx->batch_sz; b_id++) { in DMPlexLandauCreateMassMatrix()
2624 …const PetscInt moffset = LAND_MOFFSET(b_id, grid, ctx->batch_sz, ctx->num_grids, ctx->mat_offs… in DMPlexLandauCreateMassMatrix()
2645 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) PetscCall(MatDestroy(&subM[grid])); in DMPlexLandauCreateMassMatrix()
2650 ctx->M = packM; in DMPlexLandauCreateMassMatrix()
2652 PetscCall(PetscLogEventEnd(ctx->events[14], 0, 0, 0, 0)); in DMPlexLandauCreateMassMatrix()
2677 LandauCtx *ctx = (LandauCtx *)actx; in DMPlexLandauIFunction() local
2687 PetscCall(DMGetApplicationContext(pack, &ctx)); in DMPlexLandauIFunction()
2688 PetscCheck(ctx, PETSC_COMM_SELF, PETSC_ERR_PLIB, "no context"); in DMPlexLandauIFunction()
2689 if (ctx->stage) PetscCall(PetscLogStagePush(ctx->stage)); in DMPlexLandauIFunction()
2690 PetscCall(PetscLogEventBegin(ctx->events[11], 0, 0, 0, 0)); in DMPlexLandauIFunction()
2691 PetscCall(PetscLogEventBegin(ctx->events[0], 0, 0, 0, 0)); in DMPlexLandauIFunction()
2696 PetscCall(PetscObjectStateGet((PetscObject)ctx->J, &state)); in DMPlexLandauIFunction()
2697 if (state != ctx->norm_state) { in DMPlexLandauIFunction()
2698 PetscCall(MatZeroEntries(ctx->J)); in DMPlexLandauIFunction()
2699 PetscCall(LandauFormJacobian_Internal(X, ctx->J, dim, 0.0, (void *)ctx)); in DMPlexLandauIFunction()
2700 PetscCall(MatViewFromOptions(ctx->J, NULL, "-dm_landau_jacobian_view")); in DMPlexLandauIFunction()
2701 PetscCall(PetscObjectStateGet((PetscObject)ctx->J, &state)); in DMPlexLandauIFunction()
2702 ctx->norm_state = state; in DMPlexLandauIFunction()
2707 PetscCall(MatMult(ctx->J, X, F)); /* C*f */ in DMPlexLandauIFunction()
2709 if (X_t) PetscCall(MatMultAdd(ctx->M, X_t, F, F)); in DMPlexLandauIFunction()
2711 if (ctx->stage) { in DMPlexLandauIFunction()
2713 ctx->times[LANDAU_OPERATOR] += (endtime - starttime); in DMPlexLandauIFunction()
2714 ctx->times[LANDAU_JACOBIAN] += (endtime - starttime); in DMPlexLandauIFunction()
2715 ctx->times[LANDAU_MATRIX_TOTAL] += (endtime - starttime); in DMPlexLandauIFunction()
2716 ctx->times[LANDAU_JACOBIAN_COUNT] += 1; in DMPlexLandauIFunction()
2719 PetscCall(PetscLogEventEnd(ctx->events[0], 0, 0, 0, 0)); in DMPlexLandauIFunction()
2720 PetscCall(PetscLogEventEnd(ctx->events[11], 0, 0, 0, 0)); in DMPlexLandauIFunction()
2721 if (ctx->stage) PetscCall(PetscLogStagePop()); in DMPlexLandauIFunction()
2748 LandauCtx *ctx = NULL; in DMPlexLandauIJacobian() local
2758 PetscCall(DMGetApplicationContext(pack, &ctx)); in DMPlexLandauIJacobian()
2759 PetscCheck(ctx, PETSC_COMM_SELF, PETSC_ERR_PLIB, "no context"); in DMPlexLandauIJacobian()
2760 …PetscCheck(Amat == Pmat && Amat == ctx->J, ctx->comm, PETSC_ERR_PLIB, "Amat!=Pmat || Amat!=ctx->J"… in DMPlexLandauIJacobian()
2763 if (ctx->stage) PetscCall(PetscLogStagePush(ctx->stage)); in DMPlexLandauIJacobian()
2764 PetscCall(PetscLogEventBegin(ctx->events[11], 0, 0, 0, 0)); in DMPlexLandauIJacobian()
2765 PetscCall(PetscLogEventBegin(ctx->events[9], 0, 0, 0, 0)); in DMPlexLandauIJacobian()
2769 PetscCheck(shift != 0.0, ctx->comm, PETSC_ERR_PLIB, "zero shift"); in DMPlexLandauIJacobian()
2770 PetscCall(PetscObjectStateGet((PetscObject)ctx->J, &state)); in DMPlexLandauIJacobian()
2771 …PetscCheck(state == ctx->norm_state, ctx->comm, PETSC_ERR_PLIB, "wrong state, %" PetscInt64_FMT " … in DMPlexLandauIJacobian()
2772 if (!ctx->use_matrix_mass) { in DMPlexLandauIJacobian()
2773 PetscCall(LandauFormJacobian_Internal(X, ctx->J, dim, shift, (void *)ctx)); in DMPlexLandauIJacobian()
2775 PetscCall(MatAXPY(Pmat, shift, ctx->M, SAME_NONZERO_PATTERN)); in DMPlexLandauIJacobian()
2778 if (ctx->stage) { in DMPlexLandauIJacobian()
2780 ctx->times[LANDAU_OPERATOR] += (endtime - starttime); in DMPlexLandauIJacobian()
2781 ctx->times[LANDAU_MASS] += (endtime - starttime); in DMPlexLandauIJacobian()
2782 ctx->times[LANDAU_MATRIX_TOTAL] += (endtime - starttime); in DMPlexLandauIJacobian()
2785 PetscCall(PetscLogEventEnd(ctx->events[9], 0, 0, 0, 0)); in DMPlexLandauIJacobian()
2786 PetscCall(PetscLogEventEnd(ctx->events[11], 0, 0, 0, 0)); in DMPlexLandauIJacobian()
2787 if (ctx->stage) PetscCall(PetscLogStagePop()); in DMPlexLandauIJacobian()