[petsc-users] asymmetric errors at X+/Y+ boundary on periodic box mesh
Yann Jobic
yann.jobic at univ-amu.fr
Fri May 22 08:50:21 CDT 2026
Dear all,
I'm writing a small finite-volume scalar advection code on a 2D periodic
DMPlex hex mesh, using the same pipeline as src/ts/tutorials/ex11.c
(PetscFV + PetscLimiter + DMPlexTSComputeRHSFunctionFVM + TSSSP).
After only 10 time steps from a smooth periodic initial condition
phi_0(x,y) = sin(2 pi x) * sin(2 pi y), with constant velocity a =
(1,1), on a 32x32 hex periodic box, the error |phi - phi_exact| is
strongly asymmetric: it is roughly 0.18 in the cells touching x=1 and
y=1, but about a tenth of that everywhere else (including near x=0 and
y=0). The exact solution is just phi_0 translated by a*t and stays smooth.
I'm using PETSc 3.24.4 (release) and 3.25.1 (release).
I'm attaching a 170-line MWE (mwe.c, no command-line options needed) and
a screenshot of the problem. The MWE prints
t=1.000e-02 Linf=1.7897e-01
Is there an additional setup step that I'm missing for the periodic hex
case (e.g. an explicit call about face localization, or about the SF
after DMPlexConstructGhostCells)?
Thanks,
Yann
-------------- next part --------------
/* MWE: asymmetric error at X+/Y+ on a periodic hex DMPlex with PetscFV
leastsquares + vanleer + DMPlexTSComputeRHSFunctionFVM.
Build: make mwe (or compile against PETSc with -lpetsc)
Run: mpiexec -n 1 ./mwe
(writes mwe.vtu; open in ParaView, color by "error")
On PETSc 3.24.4 and 3.25.1 (release builds), the "error" field at t=0.01
is ~0.18 in cells touching x=1 and y=1, vs ~0.01 elsewhere. Velocity is
(1,1), so the exact solution at t=0.01 is just sin(2 pi (x-0.01)) *
sin(2 pi (y-0.01)). All other CLI options are passed in the code below. */
#include <petscts.h>
#include <petscdmplex.h>
#include <petscfv.h>
#include <petscds.h>
static const PetscReal A[2] = {1.0, 1.0}; /* advection velocity */
static PetscErrorCode Sine(PetscInt dim, PetscReal t, const PetscReal x[],
PetscInt Nf, PetscScalar u[], void *ctx)
{
u[0] = PetscSinReal(2*PETSC_PI*(x[0] - A[0]*t)) *
PetscSinReal(2*PETSC_PI*(x[1] - A[1]*t));
return PETSC_SUCCESS;
}
static void Riemann(PetscInt dim, PetscInt Nf, const PetscReal qp[],
const PetscReal n[], const PetscScalar uL[],
const PetscScalar uR[], PetscInt nC,
const PetscScalar constants[], PetscScalar flux[], void *ctx)
{
PetscReal wn = A[0]*n[0] + A[1]*n[1];
flux[0] = (wn > 0 ? PetscRealPart(uL[0]) : PetscRealPart(uR[0])) * wn;
}
int main(int argc, char **argv)
{
DM dm;
PetscFV fvm;
PetscDS ds;
TS ts;
Vec U;
PetscCall(PetscInitialize(&argc, &argv, NULL, NULL));
/* Mesh: periodic 32x32 hex box, FV adjacency, overlap 1. */
const char *opts =
"-dm_plex_box_faces 32,32 -dm_plex_simplex 0 "
"-dm_plex_box_bd periodic,periodic -dm_distribute_overlap 1 "
"-dm_plex_adj_cone -dm_plex_adj_closure 0";
PetscCall(PetscOptionsInsertString(NULL, opts));
PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
PetscCall(DMSetType(dm, DMPLEX));
PetscCall(DMSetFromOptions(dm));
{
DM g;
PetscCall(DMPlexConstructGhostCells(dm, NULL, NULL, &g));
PetscCall(DMDestroy(&dm));
dm = g;
}
/* FV: leastsquares + vanleer. */
PetscCall(PetscFVCreate(PETSC_COMM_WORLD, &fvm));
PetscCall(PetscFVSetType(fvm, PETSCFVLEASTSQUARES));
PetscCall(PetscFVSetNumComponents(fvm, 1));
PetscCall(PetscFVSetSpatialDimension(fvm, 2));
PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
{
PetscLimiter lim;
PetscCall(PetscFVGetLimiter(fvm, &lim));
PetscCall(PetscLimiterSetType(lim, PETSCLIMITERVANLEER));
}
PetscCall(DMAddField(dm, NULL, (PetscObject)fvm));
PetscCall(DMCreateDS(dm));
PetscCall(DMGetDS(dm, &ds));
PetscCall(PetscDSSetRiemannSolver(ds, 0, Riemann));
/* IC + TS. */
PetscCall(DMCreateGlobalVector(dm, &U));
PetscCall(PetscObjectSetName((PetscObject)U, "phi"));
{
PetscErrorCode (*f[1])(PetscInt, PetscReal, const PetscReal[],
PetscInt, PetscScalar[], void*) = {Sine};
void *fctx[1] = {NULL};
PetscCall(DMProjectFunction(dm, 0.0, f, fctx, INSERT_ALL_VALUES, U));
}
PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
PetscCall(TSSetType(ts, TSSSP));
PetscCall(TSSSPSetType(ts, TSSSPRKS3));
PetscCall(TSSetDM(ts, dm));
PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, NULL));
PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, NULL));
PetscCall(TSSetTimeStep(ts, 1e-3));
PetscCall(TSSetMaxTime(ts, 0.01));
PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
PetscCall(TSSetSolution(ts, U));
PetscCall(TSSolve(ts, U));
/* Build an "error" vector for VTK, and report Linf. */
Vec E;
PetscCall(DMCreateGlobalVector(dm, &E));
PetscCall(PetscObjectSetName((PetscObject)E, "error"));
{
Vec Uloc;
Vec cg;
DM dmCell;
PetscReal tF;
PetscInt cS, cE;
PetscSection gsec;
const PetscScalar *cgArr, *uArr;
PetscScalar *eArr;
PetscCall(TSGetTime(ts, &tF));
PetscCall(DMGetLocalVector(dm, &Uloc));
PetscCall(DMGlobalToLocalBegin(dm, U, INSERT_VALUES, Uloc));
PetscCall(DMGlobalToLocalEnd(dm, U, INSERT_VALUES, Uloc));
PetscCall(DMPlexGetGeometryFVM(dm, NULL, &cg, NULL));
PetscCall(VecGetDM(cg, &dmCell));
PetscCall(DMGetGlobalSection(dm, &gsec));
PetscCall(DMPlexGetHeightStratum(dm, 0, &cS, &cE));
PetscCall(VecGetArrayRead(cg, &cgArr));
PetscCall(VecGetArrayRead(Uloc, &uArr));
PetscCall(VecGetArray(E, &eArr));
PetscReal linf = 0.0;
for (PetscInt c = cS; c < cE; ++c) {
PetscInt gOff;
PetscCall(PetscSectionGetOffset(gsec, c, &gOff));
if (gOff < 0) continue;
PetscFVCellGeom *cgC;
const PetscScalar *uC;
PetscScalar *eC;
PetscCall(DMPlexPointLocalRead(dmCell, c, cgArr, &cgC));
PetscCall(DMPlexPointLocalRead(dm, c, uArr, &uC));
PetscCall(DMPlexPointGlobalRef(dm, c, eArr, &eC));
if (!cgC || !uC || !eC) continue;
PetscScalar uref;
Sine(2, tF, cgC->centroid, 1, &uref, NULL);
PetscReal e = PetscAbsReal(PetscRealPart(uC[0] - uref));
eC[0] = e;
if (e > linf) linf = e;
}
PetscCall(VecRestoreArray(E, &eArr));
PetscCall(VecRestoreArrayRead(Uloc, &uArr));
PetscCall(VecRestoreArrayRead(cg, &cgArr));
PetscCall(DMRestoreLocalVector(dm, &Uloc));
PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &linf, 1, MPIU_REAL, MPIU_MAX,
PETSC_COMM_WORLD));
PetscCall(PetscPrintf(PETSC_COMM_WORLD,
"t=%.3e Linf=%.4e\n", (double)tF, (double)linf));
}
/* Write the VTK so the asymmetric layer at x=1 / y=1 is visible. */
{
PetscViewer v;
PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &v));
PetscCall(PetscViewerSetType(v, PETSCVIEWERVTK));
PetscCall(PetscViewerFileSetMode(v, FILE_MODE_WRITE));
PetscCall(PetscViewerFileSetName(v, "mwe.vtu"));
PetscCall(VecView(U, v));
PetscCall(VecView(E, v));
PetscCall(PetscViewerDestroy(&v));
}
PetscCall(VecDestroy(&U));
PetscCall(VecDestroy(&E));
PetscCall(TSDestroy(&ts));
PetscCall(DMDestroy(&dm));
PetscCall(PetscFVDestroy(&fvm));
PetscCall(PetscFinalize());
return 0;
}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: t_01_phi_error.jpeg
Type: image/jpeg
Size: 17089 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20260522/0c8b257d/attachment-0001.jpeg>
More information about the petsc-users
mailing list