Skip to content

ex - add PETSc BP1+3 and BP2+4 examples #1820

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

Merged
jeremylt merged 1 commit into from
May 13, 2025
Merged

ex - add PETSc BP1+3 and BP2+4 examples #1820

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
11 changes: 8 additions & 3 deletions examples/petsc/bps.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,8 @@
//
//TESTARGS(name="BP3, tet elements") -ceed {ceed_resource} -test -problem bp3 -degree 3 -ksp_max_it_clip 50,50 -simplex
//TESTARGS(name="BP5, hex elements") -ceed {ceed_resource} -test -problem bp5 -degree 3 -ksp_max_it_clip 18,18
//TESTARGS(name="BP1+3, hex elements") -ceed {ceed_resource} -test -problem bp1_3 -degree 3 -ksp_max_it_clip 18,18
//TESTARGS(name="BP2+4, hex elements") -ceed {ceed_resource} -test -problem bp2_4 -degree 3 -ksp_max_it_clip 18,18

/// @file
/// CEED BPs example using PETSc with DMPlex
Expand Down Expand Up @@ -183,9 +185,9 @@ static PetscErrorCode RunWithDM(RunParams rp, DM dm, const char *ceed_resource)
{
PC pc;
PetscCall(KSPGetPC(ksp, &pc));
if (rp->bp_choice == CEED_BP1 || rp->bp_choice == CEED_BP2) {
if (rp->bp_choice == CEED_BP1 || rp->bp_choice == CEED_BP2 || rp->bp_choice == CEED_BP13 || rp->bp_choice == CEED_BP24) {
PetscCall(PCSetType(pc, PCJACOBI));
if (rp->simplex) {
if (rp->simplex || rp->bp_choice == CEED_BP13 || rp->bp_choice == CEED_BP24) {
PetscCall(PCJacobiSetType(pc, PC_JACOBI_DIAGONAL));
} else {
PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM));
Expand Down Expand Up @@ -255,7 +257,10 @@ static PetscErrorCode RunWithDM(RunParams rp, DM dm, const char *ceed_resource)
PetscCall(SetupErrorOperatorCtx(rp->comm, dm, ceed, ceed_data, X_loc, op_error, op_error_ctx));
PetscScalar l2_error;
PetscCall(ComputeL2Error(X, &l2_error, op_error_ctx));
PetscReal tol = 5e-2;
// Tighter tol for BP1, BP2
// Looser tol for BP3, BP4, BP5, and BP6 with extra for vector valued problems
// BP1+3 and BP2+4 follow the pattern for BP3 and BP4
PetscReal tol = rp->bp_choice < CEED_BP3 ? 5e-4 : (5e-2 + (rp->bp_choice % 2 == 1 ? 5e-3 : 0));
if (!rp->test_mode || l2_error > tol) {
PetscCall(MPI_Allreduce(&my_rt, &rt_min, 1, MPI_DOUBLE, MPI_MIN, rp->comm));
PetscCall(MPI_Allreduce(&my_rt, &rt_max, 1, MPI_DOUBLE, MPI_MAX, rp->comm));
Expand Down
2 changes: 1 addition & 1 deletion examples/petsc/bps.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,4 +17,4 @@ static const char *const mem_types[] = {"host", "device", "memType", "CEED_MEM_"
typedef enum { COARSEN_UNIFORM = 0, COARSEN_LOGARITHMIC = 1 } CoarsenType;
static const char *const coarsen_types[] = {"uniform", "logarithmic", "CoarsenType", "COARSEN", 0};

static const char *const bp_types[] = {"bp1", "bp2", "bp3", "bp4", "bp5", "bp6", "BPType", "CEED_BP", 0};
static const char *const bp_types[] = {"bp1", "bp2", "bp3", "bp4", "bp5", "bp6", "bp1_3", "bp2_4", "BPType", "CEED_BP", 0};
242 changes: 139 additions & 103 deletions examples/petsc/include/bpsproblemdata.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,9 @@

#include "../include/structs.h"
#include "../qfunctions/bps/bp1.h"
#include "../qfunctions/bps/bp13.h"
#include "../qfunctions/bps/bp2.h"
#include "../qfunctions/bps/bp24.h"
#include "../qfunctions/bps/bp3.h"
#include "../qfunctions/bps/bp4.h"
#include "../qfunctions/bps/common.h"
Expand All @@ -23,107 +25,141 @@
// BP Option Data
// -----------------------------------------------------------------------------

BPData bp_options[6] = {
[CEED_BP1] = {.num_comp_u = 1,
.num_comp_x = 3,
.topo_dim = 3,
.q_data_size = 1,
.q_extra = 1,
.setup_geo = SetupMassGeo,
.setup_rhs = SetupMassRhs,
.apply = Mass,
.error = Error,
.setup_geo_loc = SetupMassGeo_loc,
.setup_rhs_loc = SetupMassRhs_loc,
.apply_loc = Mass_loc,
.error_loc = Error_loc,
.in_mode = CEED_EVAL_INTERP,
.out_mode = CEED_EVAL_INTERP,
.q_mode = CEED_GAUSS,
.enforce_bc = PETSC_FALSE},
[CEED_BP2] = {.num_comp_u = 3,
.num_comp_x = 3,
.topo_dim = 3,
.q_data_size = 1,
.q_extra = 1,
.setup_geo = SetupMassGeo,
.setup_rhs = SetupMassRhs3,
.apply = Mass3,
.error = Error3,
.setup_geo_loc = SetupMassGeo_loc,
.setup_rhs_loc = SetupMassRhs3_loc,
.apply_loc = Mass3_loc,
.error_loc = Error3_loc,
.in_mode = CEED_EVAL_INTERP,
.out_mode = CEED_EVAL_INTERP,
.q_mode = CEED_GAUSS,
.enforce_bc = PETSC_FALSE},
[CEED_BP3] = {.num_comp_u = 1,
.num_comp_x = 3,
.topo_dim = 3,
.q_data_size = 7,
.q_extra = 1,
.setup_geo = SetupDiffGeo,
.setup_rhs = SetupDiffRhs,
.apply = Diff,
.error = Error,
.setup_geo_loc = SetupDiffGeo_loc,
.setup_rhs_loc = SetupDiffRhs_loc,
.apply_loc = Diff_loc,
.error_loc = Error_loc,
.in_mode = CEED_EVAL_GRAD,
.out_mode = CEED_EVAL_GRAD,
.q_mode = CEED_GAUSS,
.enforce_bc = PETSC_TRUE },
[CEED_BP4] = {.num_comp_u = 3,
.num_comp_x = 3,
.topo_dim = 3,
.q_data_size = 7,
.q_extra = 1,
.setup_geo = SetupDiffGeo,
.setup_rhs = SetupDiffRhs3,
.apply = Diff3,
.error = Error3,
.setup_geo_loc = SetupDiffGeo_loc,
.setup_rhs_loc = SetupDiffRhs3_loc,
.apply_loc = Diff3_loc,
.error_loc = Error3_loc,
.in_mode = CEED_EVAL_GRAD,
.out_mode = CEED_EVAL_GRAD,
.q_mode = CEED_GAUSS,
.enforce_bc = PETSC_TRUE },
[CEED_BP5] = {.num_comp_u = 1,
.num_comp_x = 3,
.topo_dim = 3,
.q_data_size = 7,
.q_extra = 0,
.setup_geo = SetupDiffGeo,
.setup_rhs = SetupDiffRhs,
.apply = Diff,
.error = Error,
.setup_geo_loc = SetupDiffGeo_loc,
.setup_rhs_loc = SetupDiffRhs_loc,
.apply_loc = Diff_loc,
.error_loc = Error_loc,
.in_mode = CEED_EVAL_GRAD,
.out_mode = CEED_EVAL_GRAD,
.q_mode = CEED_GAUSS_LOBATTO,
.enforce_bc = PETSC_TRUE },
[CEED_BP6] = {.num_comp_u = 3,
.num_comp_x = 3,
.topo_dim = 3,
.q_data_size = 7,
.q_extra = 0,
.setup_geo = SetupDiffGeo,
.setup_rhs = SetupDiffRhs3,
.apply = Diff3,
.error = Error3,
.setup_geo_loc = SetupDiffGeo_loc,
.setup_rhs_loc = SetupDiffRhs3_loc,
.apply_loc = Diff3_loc,
.error_loc = Error3_loc,
.in_mode = CEED_EVAL_GRAD,
.out_mode = CEED_EVAL_GRAD,
.q_mode = CEED_GAUSS_LOBATTO,
.enforce_bc = PETSC_TRUE }
BPData bp_options[8] = {
[CEED_BP1] = {.num_comp_u = 1,
.num_comp_x = 3,
.topo_dim = 3,
.q_data_size = 1,
.q_extra = 1,
.setup_geo = SetupMassGeo,
.setup_rhs = SetupMassRhs,
.apply = Mass,
.error = Error,
.setup_geo_loc = SetupMassGeo_loc,
.setup_rhs_loc = SetupMassRhs_loc,
.apply_loc = Mass_loc,
.error_loc = Error_loc,
.in_mode = CEED_EVAL_INTERP,
.out_mode = CEED_EVAL_INTERP,
.q_mode = CEED_GAUSS,
.enforce_bc = PETSC_FALSE},
[CEED_BP2] = {.num_comp_u = 3,
.num_comp_x = 3,
.topo_dim = 3,
.q_data_size = 1,
.q_extra = 1,
.setup_geo = SetupMassGeo,
.setup_rhs = SetupMassRhs3,
.apply = Mass3,
.error = Error3,
.setup_geo_loc = SetupMassGeo_loc,
.setup_rhs_loc = SetupMassRhs3_loc,
.apply_loc = Mass3_loc,
.error_loc = Error3_loc,
.in_mode = CEED_EVAL_INTERP,
.out_mode = CEED_EVAL_INTERP,
.q_mode = CEED_GAUSS,
.enforce_bc = PETSC_FALSE},
[CEED_BP3] = {.num_comp_u = 1,
.num_comp_x = 3,
.topo_dim = 3,
.q_data_size = 7,
.q_extra = 1,
.setup_geo = SetupDiffGeo,
.setup_rhs = SetupDiffRhs,
.apply = Diff,
.error = Error,
.setup_geo_loc = SetupDiffGeo_loc,
.setup_rhs_loc = SetupDiffRhs_loc,
.apply_loc = Diff_loc,
.error_loc = Error_loc,
.in_mode = CEED_EVAL_GRAD,
.out_mode = CEED_EVAL_GRAD,
.q_mode = CEED_GAUSS,
.enforce_bc = PETSC_TRUE },
[CEED_BP4] = {.num_comp_u = 3,
.num_comp_x = 3,
.topo_dim = 3,
.q_data_size = 7,
.q_extra = 1,
.setup_geo = SetupDiffGeo,
.setup_rhs = SetupDiffRhs3,
.apply = Diff3,
.error = Error3,
.setup_geo_loc = SetupDiffGeo_loc,
.setup_rhs_loc = SetupDiffRhs3_loc,
.apply_loc = Diff3_loc,
.error_loc = Error3_loc,
.in_mode = CEED_EVAL_GRAD,
.out_mode = CEED_EVAL_GRAD,
.q_mode = CEED_GAUSS,
.enforce_bc = PETSC_TRUE },
[CEED_BP5] = {.num_comp_u = 1,
.num_comp_x = 3,
.topo_dim = 3,
.q_data_size = 7,
.q_extra = 0,
.setup_geo = SetupDiffGeo,
.setup_rhs = SetupDiffRhs,
.apply = Diff,
.error = Error,
.setup_geo_loc = SetupDiffGeo_loc,
.setup_rhs_loc = SetupDiffRhs_loc,
.apply_loc = Diff_loc,
.error_loc = Error_loc,
.in_mode = CEED_EVAL_GRAD,
.out_mode = CEED_EVAL_GRAD,
.q_mode = CEED_GAUSS_LOBATTO,
.enforce_bc = PETSC_TRUE },
[CEED_BP6] = {.num_comp_u = 3,
.num_comp_x = 3,
.topo_dim = 3,
.q_data_size = 7,
.q_extra = 0,
.setup_geo = SetupDiffGeo,
.setup_rhs = SetupDiffRhs3,
.apply = Diff3,
.error = Error3,
.setup_geo_loc = SetupDiffGeo_loc,
.setup_rhs_loc = SetupDiffRhs3_loc,
.apply_loc = Diff3_loc,
.error_loc = Error3_loc,
.in_mode = CEED_EVAL_GRAD,
.out_mode = CEED_EVAL_GRAD,
.q_mode = CEED_GAUSS_LOBATTO,
.enforce_bc = PETSC_TRUE },
[CEED_BP13] = {.num_comp_u = 1,
.num_comp_x = 3,
.topo_dim = 3,
.q_data_size = 7,
.q_extra = 1,
.setup_geo = SetupDiffGeo,
.setup_rhs = SetupMassDiffRhs,
.apply = MassDiff,
.error = Error,
.setup_geo_loc = SetupDiffGeo_loc,
.setup_rhs_loc = SetupMassDiffRhs_loc,
.apply_loc = MassDiff_loc,
.error_loc = Error_loc,
.in_mode = CEED_EVAL_INTERP + CEED_EVAL_GRAD,
.out_mode = CEED_EVAL_INTERP + CEED_EVAL_GRAD,
.q_mode = CEED_GAUSS,
.enforce_bc = PETSC_TRUE },
[CEED_BP24] = {.num_comp_u = 3,
.num_comp_x = 3,
.topo_dim = 3,
.q_data_size = 7,
.q_extra = 1,
.setup_geo = SetupDiffGeo,
.setup_rhs = SetupMassDiffRhs3,
.apply = MassDiff3,
.error = Error3,
.setup_geo_loc = SetupDiffGeo_loc,
.setup_rhs_loc = SetupMassDiffRhs3_loc,
.apply_loc = MassDiff3_loc,
.error_loc = Error3_loc,
.in_mode = CEED_EVAL_INTERP + CEED_EVAL_GRAD,
.out_mode = CEED_EVAL_INTERP + CEED_EVAL_GRAD,
.q_mode = CEED_GAUSS,
.enforce_bc = PETSC_TRUE },
};
2 changes: 1 addition & 1 deletion examples/petsc/include/structs.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ typedef struct {
} BPData;

// BP options
typedef enum { CEED_BP1 = 0, CEED_BP2 = 1, CEED_BP3 = 2, CEED_BP4 = 3, CEED_BP5 = 4, CEED_BP6 = 5 } BPType;
typedef enum { CEED_BP1 = 0, CEED_BP2 = 1, CEED_BP3 = 2, CEED_BP4 = 3, CEED_BP5 = 4, CEED_BP6 = 5, CEED_BP13 = 6, CEED_BP24 = 7 } BPType;

// -----------------------------------------------------------------------------
// Parameter structure for running problems
Expand Down
74 changes: 74 additions & 0 deletions examples/petsc/qfunctions/bps/bp13.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,74 @@
// Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and other CEED contributors.
// All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
//
// SPDX-License-Identifier: BSD-2-Clause
//
// This file is part of CEED: http://github.com/ceed

/// @file
/// libCEED QFunctions for diffusion operator example using PETSc

#include <ceed/types.h>
#ifndef CEED_RUNNING_JIT_PASS
#include <math.h>
#endif

// -----------------------------------------------------------------------------
// This QFunction sets up the rhs and true solution for the problem
// -----------------------------------------------------------------------------
CEED_QFUNCTION(SetupMassDiffRhs)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
const CeedScalar *x = in[0], *w = in[1];
CeedScalar *true_soln = out[0], *rhs = out[1];

// Quadrature Point Loop
CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
const CeedScalar c[3] = {0, 1., 2.};
const CeedScalar k[3] = {1., 2., 3.};

true_soln[i] = sin(M_PI * (c[0] + k[0] * x[i + Q * 0])) * sin(M_PI * (c[1] + k[1] * x[i + Q * 1])) * sin(M_PI * (c[2] + k[2] * x[i + Q * 2]));

rhs[i] = w[i + Q * 0] * (M_PI * M_PI * (k[0] * k[0] + k[1] * k[1] + k[2] * k[2]) + 1.0) * true_soln[i];
} // End of Quadrature Point Loop
return 0;
}

// -----------------------------------------------------------------------------
// This QFunction applies the mass + diffusion operator for a scalar field.
//
// Inputs:
// u - Input vector at quadrature points
// ug - Input vector gradient at quadrature points
// q_data - Geometric factors
//
// Output:
// v - Output vector (test functions) at quadrature points
// vg - Output vector (test functions) gradient at quadrature points
// -----------------------------------------------------------------------------
CEED_QFUNCTION(MassDiff)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
const CeedScalar *u = in[0], *ug = in[1], *q_data = in[2];
CeedScalar *v = out[0], *vg = out[1];

// Quadrature Point Loop
CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
// Read spatial derivatives of u
const CeedScalar du[3] = {ug[i + Q * 0], ug[i + Q * 1], ug[i + Q * 2]};
// Read q_data (dXdxdXdx_T symmetric matrix)
const CeedScalar dXdxdXdx_T[3][3] = {
{q_data[i + 1 * Q], q_data[i + 2 * Q], q_data[i + 3 * Q]},
{q_data[i + 2 * Q], q_data[i + 4 * Q], q_data[i + 5 * Q]},
{q_data[i + 3 * Q], q_data[i + 5 * Q], q_data[i + 6 * Q]}
};

// Mass
v[i] = q_data[i + 0 * Q] * u[i];
// Diff
for (int j = 0; j < 3; j++) { // j = direction of vg
vg[i + j * Q] = (du[0] * dXdxdXdx_T[0][j] + du[1] * dXdxdXdx_T[1][j] + du[2] * dXdxdXdx_T[2][j]);
}
} // End of Quadrature Point Loop
return 0;
}
// -----------------------------------------------------------------------------
Loading
Loading