Skip to content

Add curved velocity test case to supg example. Also the tests were no… #181

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
stevendargaville wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Add curved velocity test case to supg example. Also the tests were no… #181

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
695dbb2
Add curved velocity test case to supg example. Also the tests were no…
stevendargaville Feb 13, 2026
92bddc3
Lumping is important for supg discretisations so may as well set it f…
stevendargaville Feb 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 12 additions & 4 deletions tests/Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,10 +234,13 @@ run_tests_no_load_serial:
@echo ""
@echo "Test AIRG with SUPG CG FEM in 2D"
./adv_diff_cg_supg -adv_diff_petscspace_degree 1 -dm_plex_simplex 0 -dm_refine 1 -pc_type air -pc_air_a_drop 1e-3 -pc_air_inverse_type power \
-ksp_type richardson -ksp_norm_type unpreconditioned -ksp_max_it 4
-ksp_type richardson -ksp_norm_type unpreconditioned -ksp_max_it 6 -pc_air_a_lump
@echo "Test AIRG with SUPG CG FEM in 2D with curved velocity"
./adv_diff_cg_supg -adv_diff_petscspace_degree 1 -dm_plex_simplex 0 -dm_refine 1 -pc_type air -pc_air_a_drop 1e-3 -pc_air_inverse_type power \
-ksp_type richardson -ksp_norm_type unpreconditioned -ksp_max_it 6 -curved_velocity -pc_air_a_lump
@echo "Test AIRG with SUPG CG FEM in 3D"
./adv_diff_cg_supg -adv_diff_petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 5,5,5 -pc_type air -pc_air_a_drop 1e-3 -pc_air_inverse_type power \
-ksp_type richardson -ksp_norm_type unpreconditioned -ksp_max_it 4
-ksp_type richardson -ksp_norm_type unpreconditioned -ksp_max_it 4 -pc_air_a_lump
#
@echo ""
@echo "Test AIRG with GMRES polynomials for 2D finite difference stencil"
Expand Down Expand Up @@ -486,10 +489,15 @@ run_tests_no_load_parallel:
@echo ""
@echo "Test AIRG with SUPG CG FEM in 2D in parallel"
$(MPIEXEC) -n 2 ./adv_diff_cg_supg -adv_diff_petscspace_degree 1 -dm_plex_simplex 0 -dm_refine 1 -pc_type air \
-ksp_type richardson -ksp_norm_type unpreconditioned -ksp_max_it 4 -pc_air_a_drop 1e-3 -pc_air_inverse_type power
-ksp_type richardson -ksp_norm_type unpreconditioned -ksp_max_it 6 -pc_air_a_drop 1e-3 -pc_air_inverse_type power \
-pc_air_a_lump
@echo "Test AIRG with SUPG CG FEM in 2D with curved velocity in parallel"
$(MPIEXEC) -n 2 ./adv_diff_cg_supg -adv_diff_petscspace_degree 1 -dm_plex_simplex 0 -dm_refine 1 -pc_type air -pc_air_a_drop 1e-3 -pc_air_inverse_type power \
-ksp_type richardson -ksp_norm_type unpreconditioned -ksp_max_it 6 -curved_velocity -pc_air_a_lump
@echo "Test AIRG with SUPG CG FEM in 3D in parallel"
$(MPIEXEC) -n 2 ./adv_diff_cg_supg -adv_diff_petscspace_degree 1 -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 5,5,5 \
-pc_type air -ksp_type richardson -ksp_norm_type unpreconditioned -ksp_max_it 4 -pc_air_a_drop 1e-3 -pc_air_inverse_type power
-pc_type air -ksp_type richardson -ksp_norm_type unpreconditioned -ksp_max_it 4 -pc_air_a_drop 1e-3 -pc_air_inverse_type power \
-pc_air_a_lump
#
@echo ""
@echo "Test AIRG with GMRES polynomials for 2D finite difference stencil in parallel"
Expand Down
70 changes: 58 additions & 12 deletions tests/adv_diff_cg_supg.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,9 +36,28 @@ static char help[] = "Solves steady advection-diffusion FEM problem with SUPG st

typedef struct {
PetscReal alpha; // Diffusion coefficient
PetscReal advection_velocity[3]; // Advection velocity, in 2D or 3D
PetscReal advection_velocity[3]; // Advection velocity, in 2D or 3D (for straight line case)
PetscBool curved_velocity; // Use curved velocity field if true
} AppCtx;

// Helper function to compute velocity at a point
// constants[0] is alpha, constants[1-3] are constant velocity, constants[4] is curved flag (-1.0 if curved)
static inline void GetVelocity(const PetscScalar constants[], const PetscReal x[], PetscReal v[])
{
if (constants[4] == -1.0) {
// Spatially-varying velocity field: top-left quadrant of a rotating circle (center at 1,0)
// u(x,y) = y, v(x,y) = 1-x
v[0] = x[1]; // u(x,y)
v[1] = 1.0 - x[0]; // v(x,y)
v[2] = 0.0; // w (unused in 2D)
} else {
// Constant velocity field
v[0] = constants[1];
v[1] = constants[2];
v[2] = constants[3];
}
}

// Helper function to compute the SUPG stabilization parameter tau
static inline PetscErrorCode ComputeSUPGStabilization(PetscInt dim, PetscReal h, PetscReal alpha, const PetscReal v[], PetscReal *tau)
{
Expand Down Expand Up @@ -92,9 +111,8 @@ static PetscErrorCode neumann_bc_zero_flux(PetscInt dim, PetscReal time, const P
// The f0 term in the weak form integral
static void advection_diffusion_f0(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
{
// constants[1] is u, constants[2] is v, constants[3] is w
// Pointer to the first component of the velocity in the constants array
const PetscReal *advection_velocity = &constants[1];
PetscReal advection_velocity[3];
GetVelocity(constants, x, advection_velocity);

PetscReal adv_dot_grad_u = 0.0;
PetscReal volumetric_source = 0.0; // Assuming f = 0 for now
Expand All @@ -112,8 +130,12 @@ static void advection_diffusion_f0(PetscInt dim, PetscInt Nf, PetscInt NfAux, co
static void advection_diffusion_f1_supg(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
{
const PetscReal alpha = constants[0];
const PetscReal *v = &constants[1];
const PetscReal h = a[0]; // Cell size from auxiliary field

// Get velocity field
PetscReal v[3];
GetVelocity(constants, x, v);

PetscReal tau, v_dot_grad_u = 0.0;
PetscInt d;

Expand All @@ -129,9 +151,8 @@ static void advection_diffusion_f1_supg(PetscInt dim, PetscInt Nf, PetscInt NfAu
// The Jacobian for advection term
static void g1_jacobian_advection_term(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
{
// constants[1] is u, constants[2] is v, constants[3] is w
// Pointer to the first component of the velocity
const PetscReal *advection_velocity = &constants[1];
PetscReal advection_velocity[3];
GetVelocity(constants, x, advection_velocity);

PetscInt d;
for (d = 0; d < dim; ++d) {
Expand All @@ -143,8 +164,12 @@ static void g1_jacobian_advection_term(PetscInt dim, PetscInt Nf, PetscInt NfAux
static void g3_jacobian_diffusion_term_supg(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
{
const PetscReal alpha = constants[0];
const PetscReal *v = &constants[1];
const PetscReal h = a[0]; // Cell size from auxiliary field

// Get velocity field
PetscReal v[3];
GetVelocity(constants, x, v);

PetscReal tau;
PetscInt d, c;

Expand All @@ -166,6 +191,10 @@ static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
options->alpha = 0.0;
PetscCall(PetscOptionsGetReal(NULL, NULL, "-alpha", &options->alpha, NULL));

// Curved velocity field option
options->curved_velocity = PETSC_FALSE;
PetscCall(PetscOptionsGetBool(NULL, NULL, "-curved_velocity", &options->curved_velocity, NULL));

// Initialize advection to zero
options->advection_velocity[0] = 0.0; // u
options->advection_velocity[1] = 0.0; // v
Expand Down Expand Up @@ -276,13 +305,14 @@ static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *options)

/* Setup constants that get passed into the FEM functions*/
{
PetscScalar constants[4];
PetscScalar constants[5];

constants[0] = options->alpha;
constants[1] = options->advection_velocity[0];
constants[2] = options->advection_velocity[1];
constants[3] = options->advection_velocity[2];
PetscCall(PetscDSSetConstants(ds, 4, constants));
constants[4] = options->curved_velocity ? -1.0 : 0.0; // -1.0 indicates curved velocity
PetscCall(PetscDSSetConstants(ds, 5, constants));
}

PetscFunctionReturn(PETSC_SUCCESS);
Expand Down Expand Up @@ -385,6 +415,8 @@ int main(int argc, char **argv)
Vec u; /* Solutions */
AppCtx options; /* options-defined work context */
PetscLogStage setup, gpu_copy;
KSPConvergedReason reason;
KSP ksp;

PetscFunctionBeginUser;
PetscCall(PetscInitialize(&argc, &argv, NULL, help));
Expand Down Expand Up @@ -427,14 +459,28 @@ int main(int argc, char **argv)
PetscCall(SNESSolve(snes, NULL, u));
PetscCall(PetscLogStagePop());

// Get the iteration count
PetscCall(SNESGetKSP(snes, &ksp));
PetscCall(KSPGetConvergedReason(ksp,&reason));
if (reason < 0)
{
return 1;
}

// Solve
// We set x to 1 rather than random as the vecrandom doesn't yet have a
// gpu implementation and we don't want a copy occuring back to the cpu
if (second_solve)
{
PetscCall(VecSet(u, 1.0));
PetscCall(SNESSolve(snes, NULL, u));
}
}

PetscCall(KSPGetConvergedReason(ksp,&reason));
if (reason < 0)
{
return 1;
}

PetscCall(SNESGetSolution(snes, &u));
/* Cleanup */
Expand Down
Loading