Lines Matching refs:grid
37 for (PetscInt grid = 0; grid < maps[0].numgrids; grid++) { in LandauGPUMapsDestroy() local
38 PetscCall(PetscFree(maps[grid].c_maps)); in LandauGPUMapsDestroy()
39 PetscCall(PetscFree(maps[grid].gIdx)); in LandauGPUMapsDestroy()
116 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) PetscCall(DMCreateMatrix(ctx->plex[grid], &… in LandauFormJacobian_Internal() local
127 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauFormJacobian_Internal() local
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()
131 numCells[grid] = cEnd - cStart; // grids can have different topology in LandauFormJacobian_Internal()
148 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauFormJacobian_Internal() local
149 PetscCall(DMGetLocalSection(ctx->plex[grid], §ion[grid])); in LandauFormJacobian_Internal()
150 PetscCall(DMGetGlobalSection(ctx->plex[grid], &globsection[grid])); in LandauFormJacobian_Internal()
151 PetscCall(PetscSectionGetNumFields(section[grid], &Nf[grid])); in LandauFormJacobian_Internal()
155 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) cellClosure_sz += Nb * Nf[grid] * numCells[… in LandauFormJacobian_Internal() local
163 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauFormJacobian_Internal() local
164 …Vec locX = locXArray[LAND_PACK_IDX(b_id, grid)], globX = globXArray[LAND_PACK_IDX(b_id, 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()
173 …PetscCall(PetscMemcpy(cellClosure_it, coef, Nb * Nf[grid] * sizeof(*cellClosure_it))); /* change i… in LandauFormJacobian_Internal()
174 … PetscCall(DMPlexVecRestoreClosure(ctx->plex[grid], section[grid], locX2, ei, NULL, &coef)); in LandauFormJacobian_Internal()
175 cellClosure_it += Nb * Nf[grid]; in LandauFormJacobian_Internal()
216 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauFormJacobian_Internal() local
217 PetscCall(DMGetLocalSection(ctx->plex[grid], §ion[grid])); in LandauFormJacobian_Internal()
218 PetscCall(DMGetGlobalSection(ctx->plex[grid], &globsection[grid])); in LandauFormJacobian_Internal()
219 PetscCall(PetscSectionGetNumFields(section[grid], &Nf[grid])); in LandauFormJacobian_Internal()
225 for (PetscInt grid = 0; grid < num_grids; grid++) { in LandauFormJacobian_Internal() local
226 PetscInt nfloc = ctx->species_offset[grid + 1] - ctx->species_offset[grid]; in LandauFormJacobian_Internal()
227 elem_offset[grid + 1] = elem_offset[grid] + numCells[grid]; in LandauFormJacobian_Internal()
228 ip_offset[grid + 1] = ip_offset[grid] + numCells[grid] * Nq; in LandauFormJacobian_Internal()
229 ipf_offset[grid + 1] = ipf_offset[grid] + Nq * nfloc * numCells[grid]; in LandauFormJacobian_Internal()
246 PetscInt grid = 0; in LandauFormJacobian_Internal() local
247 while (b_elem_idx >= elem_offset[grid + 1]) grid++; // yuck search for grid in LandauFormJacobian_Internal()
249 …_nip = numCells[grid] * Nq, loc_Nf = ctx->species_offset[grid + 1] - ctx->species_offset[grid], lo… in LandauFormJacobian_Internal()
250 …const PetscInt moffset = LAND_MOFFSET(b_id, grid, ctx->batch_sz, ctx->num_grids, ctx->mat_offset);… in LandauFormJacobian_Internal()
252 …PetscReal *invJe = &invJ_a[(ip_offset[grid] + loc_elem * Nq) * dim * dim]; // ingJ is static d… in LandauFormJacobian_Internal()
255 …coef = &cellClosure[b_id * IPf_sz_glb + ipf_offset[grid] + loc_elem * Nb * loc_Nf]; // this is con… in LandauFormJacobian_Internal()
259 LandauIdx *const Idxs = &maps[grid].gIdx[loc_elem][f][0]; in LandauFormJacobian_Internal()
267 for (q = 0; q < maps[grid].num_face; q++) { in LandauFormJacobian_Internal()
268 PetscInt id = maps[grid].c_maps[idx][q].gid; in LandauFormJacobian_Internal()
269 PetscScalar scale = maps[grid].c_maps[idx][q].scale; in LandauFormJacobian_Internal()
284 … const PetscInt idx = b_id * IPf_sz_glb + ipf_offset[grid] + f * loc_nip + loc_elem * Nq + qi; in LandauFormJacobian_Internal()
316 PetscInt grid = 0; in LandauFormJacobian_Internal() local
321 while (glb_elem_idx >= elem_offset[grid + 1]) grid++; in LandauFormJacobian_Internal()
323 …cInt loc_Nf = ctx->species_offset[grid + 1] - ctx->species_offset[grid], loc_elem = glb_elem_idx… in LandauFormJacobian_Internal()
324 …const PetscInt moffset = LAND_MOFFSET(b_id, grid, ctx->batch_sz, ctx->num_grids, ctx->mat_offset… in LandauFormJacobian_Internal()
326 const PetscReal *invJe = &invJ_a[(ip_offset[grid] + loc_elem * Nq) * dim * dim]; in LandauFormJacobian_Internal()
335 const PetscInt jpidx_glb = ip_offset[grid] + qj + loc_elem * Nq; in LandauFormJacobian_Internal()
369 … temp1[0] += dudx[idx] * nu_beta[f + f_off] * invMass[f + f_off] * (*lambdas)[grid][grid_r]; in LandauFormJacobian_Internal()
370 … temp1[1] += dudy[idx] * nu_beta[f + f_off] * invMass[f + f_off] * (*lambdas)[grid][grid_r]; in LandauFormJacobian_Internal()
372 … temp1[2] += dudz[idx] * nu_beta[f + f_off] * invMass[f + f_off] * (*lambdas)[grid][grid_r]; in LandauFormJacobian_Internal()
374 temp2 += ff[idx] * nu_beta[f + f_off] * (*lambdas)[grid][grid_r]; 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()
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()
487 LandauIdx *const Idxs = &maps[grid].gIdx[loc_elem][fieldA][0]; in LandauFormJacobian_Internal()
496 for (q = 0, nr = 0; q < maps[grid].num_face; q++, nr++) { in LandauFormJacobian_Internal()
497 if (maps[grid].c_maps[idx][q].gid < 0) break; in LandauFormJacobian_Internal()
498 rows0[q] = maps[grid].c_maps[idx][q].gid; in LandauFormJacobian_Internal()
499 row_scale[q] = maps[grid].c_maps[idx][q].scale; in LandauFormJacobian_Internal()
510 nc = maps[grid].num_face; in LandauFormJacobian_Internal()
511 for (q = 0, nc = 0; q < maps[grid].num_face; q++, nc++) { in LandauFormJacobian_Internal()
512 if (maps[grid].c_maps[idx][q].gid < 0) break; in LandauFormJacobian_Internal()
513 cols0[q] = maps[grid].c_maps[idx][q].gid; in LandauFormJacobian_Internal()
514 col_scale[q] = maps[grid].c_maps[idx][q].scale; in LandauFormJacobian_Internal()
554 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauFormJacobian_Internal() local
555 …const PetscInt moffset = LAND_MOFFSET(b_id, grid, ctx->batch_sz, ctx->num_grids, ctx->mat_offs… in LandauFormJacobian_Internal()
559 Mat B = subJ[LAND_PACK_IDX(b_id, grid)]; in LandauFormJacobian_Internal()
594 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauDMCreateVMeshes() local
595 PetscReal par_radius = ctx->radius_par[grid], perp_radius = ctx->radius_perp[grid]; in LandauDMCreateVMeshes()
604 …, ctx->cells0, lo, hi, periodicity, PETSC_TRUE, 0, PETSC_TRUE, &ctx->plex[grid])); // TODO: make c… 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()
616 str[21] += grid; in LandauDMCreateVMeshes()
617 …scCall(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()
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 …self, 2, numCells, numVerts, cell_size, ctx->interpolate, pcell, 2, flatCoords, &ctx->plex[grid])); 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()
684 …PetscReal rad = ctx->radius[grid], inner_rad = rad * ctx->sphere_inner_radius_90degree[grid],… in LandauDMCreateVMeshes()
716 …_size, ctx->interpolate, (const PetscInt *)cells, 3, (const PetscReal *)coords, &ctx->plex[grid])); 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()
733 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauDMCreateVMeshes() local
737 PetscCall(DMConvert(ctx->plex[grid], convType, &dmforest)); in LandauDMCreateVMeshes()
742 PetscCall(DMDestroy(&ctx->plex[grid])); in LandauDMCreateVMeshes()
743 ctx->plex[grid] = dmforest; // Forest for adaptivity in LandauDMCreateVMeshes()
753 static PetscErrorCode SetupDS(DM pack, PetscInt dim, PetscInt grid, const char prefix[], LandauCtx … in SetupDS() argument
760 for (ii = ctx->species_offset[grid], i0 = 0; ii < ctx->species_offset[grid + 1]; ii++, i0++) { 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()
849 …DM dm, Vec X, PetscReal time, PetscReal temps[], PetscReal ns[], PetscInt grid, PetscInt b_id, Pet… in DMPlexLandauAddMaxwellians() argument
859 …for (PetscInt ii = ctx->species_offset[grid], i0 = 0; ii < ctx->species_offset[grid + 1]; ii++, i0… in DMPlexLandauAddMaxwellians()
892 static PetscErrorCode LandauSetInitialCondition(DM dm, Vec X, PetscInt grid, PetscInt b_id, PetscIn… in LandauSetInitialCondition() argument
899 …PetscCall(DMPlexLandauAddMaxwellians(dm, X, 0.0, ctx->thermal_temps, ctx->n, grid, b_id, n_batch, … in LandauSetInitialCondition()
907 static PetscErrorCode adaptToleranceFEM(PetscFE fem, Vec sol, PetscInt type, PetscInt grid, LandauC… in adaptToleranceFEM() argument
917 forest = ctx->plex[grid]; in adaptToleranceFEM()
933 PetscCall(PetscInfo(sol, "%" PetscInt_FMT ") Refine phase: %s\n", grid, s_refine_names[type])); in adaptToleranceFEM()
951 … r=%e, cell %" PetscInt_FMT ", qp %" PetscInt_FMT "/%" PetscInt_FMT "\n", grid, (double)r, c, qj +… in adaptToleranceFEM()
957 … cell %" PetscInt_FMT ", qp %" PetscInt_FMT "/%" PetscInt_FMT ", d=%e\n", grid, (double)r, c, qj +… in adaptToleranceFEM()
963 … PetscInt_FMT " origin cells %" PetscInt_FMT ",%" PetscInt_FMT " r=%g\n", grid, nr, rCellIdx[0], r… in adaptToleranceFEM()
984 …PetscCall(PetscInfo(sol, "\t%" PetscInt_FMT ") (debug) found %s cells\n", grid, s_refine_names[typ… in adaptToleranceFEM()
987 …PetscCall(PetscInfo(sol, "\t%" PetscInt_FMT ") (debug) found %s cells\n", grid, s_refine_names[typ… in adaptToleranceFEM()
997 …PetscCall(PetscInfo(sol, "\t%" PetscInt_FMT ") Refined %" PetscInt_FMT " cells\n", grid, nrefined)… in adaptToleranceFEM()
1007 …") %" PetscInt_FMT " cells, %" PetscInt_FMT " total quadrature points\n", grid, cEnd - cStart, Nq … 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], (gr… 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()
1097 for (PetscInt grid = 0; grid < LANDAU_MAX_GRIDS; grid++) { in ProcessOptions() local
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()
1104 num_species_grid[grid] = 0; in ProcessOptions()
1105 ctx->plex[grid] = NULL; /* cache as expensive to Convert */ in ProcessOptions()
1126 for (PetscInt grid = 0; grid < LANDAU_NUM_TIMERS; grid++) ctx->times[grid] = 0; in ProcessOptions() local
1224 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in ProcessOptions() local
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()
1231 for (PetscInt grid = 0; grid < LANDAU_MAX_GRIDS; grid++) { in ProcessOptions() local
1232 …for (PetscInt gridj = 0; gridj < LANDAU_MAX_GRIDS; gridj++) ctx->lambdas[gridj][grid] = lnLam; /* … in ProcessOptions()
1257 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in ProcessOptions() local
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 …ummy, "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()
1290 …for (PetscInt grid = nt; grid < ctx->num_grids; grid++) ctx->sphere_inner_radius_90degree[grid] = … in ProcessOptions() local
1293 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) ctx->sphere_inner_radius_90degree[grid] = 0… in ProcessOptions() local
1296 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) ctx->sphere_inner_radius_90degree[grid] = 0… in ProcessOptions() local
1298 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) ctx->sphere_inner_radius_90degree[grid] = 0… in ProcessOptions() local
1305 …for (PetscInt grid = nt; grid < ctx->num_grids; grid++) ctx->sphere_inner_radius_45degree[grid] = … in ProcessOptions() local
1308 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) ctx->sphere_inner_radius_45degree[grid] = 0… in ProcessOptions() local
1310 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) ctx->sphere_inner_radius_45degree[grid] = 0… in ProcessOptions() local
1315 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in ProcessOptions() local
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()
1416 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in CreateStaticData() local
1417 for (PetscInt ii = ctx->species_offset[grid]; ii < ctx->species_offset[grid + 1]; ii++) { in CreateStaticData()
1425 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in CreateStaticData() local
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() local
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() local
1436 PetscInt iii = ctx->species_offset[grid]; in CreateStaticData()
1437 for (PetscInt ii = iii; ii < ctx->species_offset[grid + 1]; ii++) { in CreateStaticData()
1440 …tf(PETSC_COMM_WORLD, " %14.9e", (double)(nu_alpha[ii] * nu_beta[jj] * ctx->lambdas[grid][gridj]))); in CreateStaticData()
1446 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in CreateStaticData() local
1447 PetscInt iii = ctx->species_offset[grid]; in CreateStaticData()
1448 for (PetscInt ii = iii; ii < ctx->species_offset[grid + 1]; ii++) { in CreateStaticData()
1451 …+ 1]; jj++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %14.9e", (double)ctx->lambdas[grid][gridj])); in CreateStaticData()
1465 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in CreateStaticData() local
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()
1469 numCells[grid] = cEnd - cStart; // grids can have different topology in CreateStaticData()
1470 PetscCall(DMGetLocalSection(ctx->plex[grid], §ion[grid])); in CreateStaticData()
1471 PetscCall(DMGetGlobalSection(ctx->plex[grid], &globsection[grid])); in CreateStaticData()
1472 PetscCall(PetscSectionGetNumFields(section[grid], &Nf[grid])); in CreateStaticData()
1473 ncellsTot += numCells[grid]; in CreateStaticData()
1501 for (PetscInt grid = 0, glb_elem_idx = 0; grid < ctx->num_grids; grid++) { in CreateStaticData() local
1502 PetscInt cStart, cEnd, Nfloc = Nf[grid], totDim = Nfloc * Nb; in CreateStaticData()
1503 if (grid_batch_is_inv[grid]) PetscCall(ISGetIndices(grid_batch_is_inv[grid], &plex_batch)); in CreateStaticData()
1505 PetscCall(DMPlexGetHeightStratum(ctx->plex[grid], 0, &cStart, &cEnd)); in CreateStaticData()
1507 maps[grid].d_self = NULL; in CreateStaticData()
1508 maps[grid].num_elements = numCells[grid]; in CreateStaticData()
1509 maps[grid].num_face = (PetscInt)(pow(Nq, 1. / ((double)dim)) + .001); // Q in CreateStaticData()
1510 … maps[grid].num_face = (PetscInt)(pow(maps[grid].num_face, (double)(dim - 1)) + .001); // Q^2 in CreateStaticData()
1511 maps[grid].num_reduced = 0; in CreateStaticData()
1512 maps[grid].deviceType = ctx->deviceType; in CreateStaticData()
1513 maps[grid].numgrids = ctx->num_grids; in CreateStaticData()
1515 PetscCall(PetscMalloc(maps[grid].num_elements * sizeof(*maps[grid].gIdx), &maps[grid].gIdx)); in CreateStaticData()
1518 for (PetscInt fieldA = 0; fieldA < Nf[grid]; fieldA++) { in CreateStaticData()
1525 …PetscCall(DMPlexGetClosureIndices(ctx->plex[grid], section[grid], globsection[grid], ej, PETSC_TRU… in CreateStaticData()
1528 … for (PetscInt q = 0; q < numindices; ++q) maps[grid].gIdx[eidx][fieldA][q] = indices[q]; in CreateStaticData()
1536 maps[grid].gIdx[eidx][fieldA][q] = plex_batch[indices[f]]; in CreateStaticData()
1538 maps[grid].gIdx[eidx][fieldA][q] = indices[f]; in CreateStaticData()
1545 …maps[grid].gIdx[eidx][fieldA][q] = -maps[grid].num_reduced - 1; // store (-)index: id = -(idx+1): … in CreateStaticData()
1549 …for (ii = 0, pointMaps[maps[grid].num_reduced][jj].scale = 0; ii < maps[grid].num_face; ii++) { //… in CreateStaticData()
1551 … pointMaps[maps[grid].num_reduced][jj].scale += PetscRealPart(elMat[f * numindices + ff + ii]); in CreateStaticData()
1554 sum += pointMaps[maps[grid].num_reduced][jj].scale; // diagnostic in CreateStaticData()
1556 …if (pointMaps[maps[grid].num_reduced][jj].scale == 0) pointMaps[maps[grid].num_reduced][jj].gid = … in CreateStaticData()
1559 pointMaps[maps[grid].num_reduced][jj].gid = plex_batch[indices[f]]; in CreateStaticData()
1561 pointMaps[maps[grid].num_reduced][jj].gid = indices[f]; in CreateStaticData()
1565 … } while (++jj < maps[grid].num_face && ++f < numindices); // jj is incremented if we hit the end in CreateStaticData()
1566 while (jj < maps[grid].num_face) { in CreateStaticData()
1567 pointMaps[maps[grid].num_reduced][jj].scale = 0; in CreateStaticData()
1568 pointMaps[maps[grid].num_reduced][jj].gid = -1; in CreateStaticData()
1575 …face=%" PetscInt_FMT ")\n", eidx, q, fieldA, (double)sum, LANDAU_MAX_Q_FACE, maps[grid].num_face)); in CreateStaticData()
1582 maps[grid].num_reduced++; in CreateStaticData()
1583 …ck(maps[grid].num_reduced < MAP_BF_SIZE, PETSC_COMM_SELF, PETSC_ERR_PLIB, "maps[grid].num_reduced … 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()
1603 … PetscCall(PetscMalloc(maps[grid].num_reduced * sizeof(*maps[grid].c_maps), &maps[grid].c_maps)); in CreateStaticData()
1604 for (PetscInt ej = 0; ej < maps[grid].num_reduced; ++ej) { in CreateStaticData()
1605 for (PetscInt q = 0; q < maps[grid].num_face; ++q) { in CreateStaticData()
1606 maps[grid].c_maps[ej][q].scale = pointMaps[ej][q].scale; in CreateStaticData()
1607 maps[grid].c_maps[ej][q].gid = pointMaps[ej][q].gid; in CreateStaticData()
1611 …= LANDAU_KOKKOS) PetscCall(LandauKokkosCreateMatMaps(maps, pointMaps, Nf, grid)); // implies Kokko… in CreateStaticData()
1614 PetscCall(ISRestoreIndices(grid_batch_is_inv[grid], &plex_batch)); in CreateStaticData()
1615 PetscCall(ISDestroy(&grid_batch_is_inv[grid])); // we are done with this in CreateStaticData()
1625 for (PetscInt grid = 0, glb_elem_idx = 0; grid < ctx->num_grids; grid++) { in CreateStaticData() local
1626 for (PetscInt ej = 0; ej < numCells[grid]; ++ej, glb_elem_idx++) { in CreateStaticData()
1628 …const LandauIdx *const Idxs = &maps[grid].gIdx[ej][0][0]; // just use field-0 maps, Th… in CreateStaticData()
1638 for (PetscInt q = 0; q < maps[grid].num_face; q++) { in CreateStaticData()
1639 if (maps[grid].c_maps[idx][q].gid < 0) break; in CreateStaticData()
1651 for (PetscInt grid = 0, glb_elem_idx = 0; grid < ctx->num_grids; grid++) { in CreateStaticData() local
1652 … const PetscInt moffset = LAND_MOFFSET(b_id, grid, ctx->batch_sz, ctx->num_grids, ctx->mat_offset); in CreateStaticData()
1653 for (PetscInt ej = 0; ej < numCells[grid]; ++ej, glb_elem_idx++) { in CreateStaticData()
1656 for (PetscInt fieldA = 0; fieldA < Nf[grid]; fieldA++) { in CreateStaticData()
1657 const LandauIdx *const Idxs = &maps[grid].gIdx[ej][fieldA][0]; in CreateStaticData()
1664 for (PetscInt q = 0; q < nr; q++) rows[q] = maps[grid].c_maps[idx][q].gid; in CreateStaticData()
1671 for (PetscInt q = 0; q < nc; q++) cols[q] = maps[grid].c_maps[idx][q].gid; in CreateStaticData()
1700 PetscInt outer_ipidx, outer_ej, grid, nip_glb = 0; in CreateStaticData() local
1704 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) nip_glb += Nq * numCells[grid]; in CreateStaticData() local
1711 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) ncells += numCells[grid]; in CreateStaticData() local
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()
1735 PetscCall(DMClone(ctx->plex[grid], &dmEnergy)); in CreateStaticData()
1753 for (ej = 0; ej < numCells[grid]; ++ej, ++outer_ej) { in CreateStaticData()
1757 …PetscCall(DMPlexComputeCellGeometryFEM(ctx->plex[grid], ej + cStart, quad, vj, Jdummy, invJe, detJ… in CreateStaticData()
1840 for (PetscInt grid = 0; grid < LANDAU_MAX_GRIDS; grid++) { in CreateStaticData() local
1842 …nt gridj = 0; gridj < LANDAU_MAX_GRIDS; gridj++) (*lambdas)[grid][gridj] = ctx->lambdas[grid][grid… in CreateStaticData()
1885 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauCreateJacobianMatrix() local
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()
1907 subM[grid] = gMat; in LandauCreateJacobianMatrix()
1912 PetscCall(ISInvertPermutation(isrow, PETSC_DECIDE, &grid_batch_is_inv[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()
1930 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in LandauCreateJacobianMatrix() local
1931 Mat B = subM[grid]; 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()
1956 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) PetscCall(MatDestroy(&subM[grid])); in LandauCreateJacobianMatrix() local
2068 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauCreateVelocitySpace() local
2070 PetscCall(SetupDS(ctx->plex[grid], dim, grid, prefix, ctx)); in DMPlexLandauCreateVelocitySpace()
2072 PetscCall(DMCreateGlobalVector(ctx->plex[grid], &Xsub[grid])); in DMPlexLandauCreateVelocitySpace()
2073 PetscCall(PetscObjectSetName((PetscObject)Xsub[grid], "u_orig")); in DMPlexLandauCreateVelocitySpace()
2075 PetscCall(LandauSetInitialCondition(ctx->plex[grid], Xsub[grid], grid, 0, 1, ctx)); 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()
2085 …scCall(LandauSphereMesh(ctx->plex[grid], ctx->radius[grid] * ctx->sphere_inner_radius_90degree[gri… in DMPlexLandauCreateVelocitySpace()
2086 PetscCall(LandauSetInitialCondition(ctx->plex[grid], Xsub[grid], grid, 0, 1, ctx)); in DMPlexLandauCreateVelocitySpace()
2088 if (grid == 0) { in DMPlexLandauCreateVelocitySpace()
2089 PetscCall(DMViewFromOptions(ctx->plex[grid], NULL, "-dm_landau_amr_dm_view")); in DMPlexLandauCreateVelocitySpace()
2090 PetscCall(VecSetOptionsPrefix(Xsub[grid], prefix)); in DMPlexLandauCreateVelocitySpace()
2091 PetscCall(VecViewFromOptions(Xsub[grid], NULL, "-dm_landau_amr_vec_view")); in DMPlexLandauCreateVelocitySpace()
2094 PetscCall(DMCompositeAddDM(*pack, ctx->plex[grid])); in DMPlexLandauCreateVelocitySpace()
2097 PetscCall(DMCompositeAddDM(*pack, ctx->plex[grid])); in DMPlexLandauCreateVelocitySpace()
2100 PetscCall(DMSetApplicationContext(ctx->plex[grid], ctx)); in DMPlexLandauCreateVelocitySpace()
2105 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) PetscCall(DMCompositeAddDM(*pack, ctx->plex… in DMPlexLandauCreateVelocitySpace() local
2110 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauCreateVelocitySpace() local
2112 PetscCall(VecGetLocalSize(Xsub[grid], &n)); in DMPlexLandauCreateVelocitySpace()
2113 ctx->mat_offset[grid + 1] = ctx->mat_offset[grid] + n; in DMPlexLandauCreateVelocitySpace()
2125 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauCreateVelocitySpace() local
2127 PetscCall(VecGetLocalSize(Xsub[grid], &n)); 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()
2132 PetscCall(VecGetArrayRead(Xsub[grid], &values)); // Drop whole grid in Plex ordering in DMPlexLandauCreateVelocitySpace()
2134 PetscCall(VecRestoreArrayRead(Xsub[grid], &values)); in DMPlexLandauCreateVelocitySpace()
2138 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) PetscCall(VecDestroy(&Xsub[grid])); in DMPlexLandauCreateVelocitySpace() local
2223 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauAccess() local
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()
2236 … const PetscInt moffset = LAND_MOFFSET(b_id, grid, ctx->batch_sz, ctx->num_grids, ctx->mat_offset); in DMPlexLandauAccess()
2239 PetscCall(func(vdm, vec, i0, grid, b_id, user_ctx)); in DMPlexLandauAccess()
2310 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) PetscCall(DMDestroy(&ctx->plex[grid])); in DMPlexLandauDestroyVelocitySpace() local
2414 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauPrintNorms() local
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()
2422 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2425 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2428 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2436 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2441 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2445 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2449 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2453 if (ii == 0 && grid == 0) { // do all at once in DMPlexLandauPrintNorms()
2466 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { // yes a grid loop in a grid loop to prin… in DMPlexLandauPrintNorms() local
2467 Vec v1 = globGammaArray[LAND_PACK_IDX(ctx->batch_view_idx, grid)]; 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() local
2477 PetscInt Nf = ctx->species_offset[grid + 1] - ctx->species_offset[grid], N, bs; in DMPlexLandauPrintNorms()
2478 …fArray[LAND_PACK_IDX(ctx->batch_view_idx, grid)], Gsub = globGammaArray[LAND_PACK_IDX(ctx->batch_v… in DMPlexLandauPrintNorms()
2492 for (PetscInt i = 0, ix = ctx->species_offset[grid]; i < Nf; i++, ix++) { in DMPlexLandauPrintNorms()
2511 PetscCall(DMPlexComputeIntegralFEM(ctx->plex[grid], Xloc, tt, ctx)); in DMPlexLandauPrintNorms()
2574 for (PetscInt grid = 0, ix = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauCreateMassMatrix() local
2575 PetscCall(DMClone(ctx->plex[grid], &massDM[grid])); in DMPlexLandauCreateMassMatrix()
2576 PetscCall(DMCopyFields(ctx->plex[grid], PETSC_DETERMINE, PETSC_DETERMINE, massDM[grid])); in DMPlexLandauCreateMassMatrix()
2577 PetscCall(DMCreateDS(massDM[grid])); in DMPlexLandauCreateMassMatrix()
2578 PetscCall(DMGetDS(massDM[grid], &prob)); in DMPlexLandauCreateMassMatrix()
2579 for (ix = 0, ii = ctx->species_offset[grid]; ii < ctx->species_offset[grid + 1]; ii++, ix++) { in DMPlexLandauCreateMassMatrix()
2584 PetscCall(DMCompositeAddDM(mass_pack, massDM[grid])); in DMPlexLandauCreateMassMatrix()
2587 PetscCall(DMCompositeAddDM(mass_pack, massDM[grid])); in DMPlexLandauCreateMassMatrix()
2590 PetscCall(DMCreateMatrix(massDM[grid], &subM[grid])); in DMPlexLandauCreateMassMatrix()
2595 …for (PetscInt grid = 0; grid < ctx->num_grids; grid++) PetscCall(DMCompositeAddDM(mass_pack, massD… in DMPlexLandauCreateMassMatrix() local
2605 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauCreateMassMatrix() local
2607 DM plex = massDM[grid]; in DMPlexLandauCreateMassMatrix()
2610 PetscCall(DMPlexSNESComputeJacobianFEM(plex, locX, subM[grid], subM[grid], ctx)); in DMPlexLandauCreateMassMatrix()
2612 PetscCall(DMDestroy(&massDM[grid])); in DMPlexLandauCreateMassMatrix()
2618 for (PetscInt grid = 0; grid < ctx->num_grids; grid++) { in DMPlexLandauCreateMassMatrix() local
2619 Mat B = subM[grid]; 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() local