From 07b3e9edc999cf7ed509768180b31cba99d9198a Mon Sep 17 00:00:00 2001 From: valeriabarra Date: Mon, 6 Jul 2020 18:35:30 -0600 Subject: [PATCH] SWE: Add advection test case and skeleton for geostrophic test --- .../shallow-water/qfunctions/advection.h | 433 ++++++++++++++++++ .../{shallowwater.h => geostrophic.h} | 75 ++- examples/fluids/shallow-water/shallowwater.c | 75 ++- examples/fluids/shallow-water/src/setup.c | 73 ++- examples/fluids/shallow-water/sw_headers.h | 72 ++- 5 files changed, 637 insertions(+), 91 deletions(-) create mode 100644 examples/fluids/shallow-water/qfunctions/advection.h rename examples/fluids/shallow-water/qfunctions/{shallowwater.h => geostrophic.h} (87%) diff --git a/examples/fluids/shallow-water/qfunctions/advection.h b/examples/fluids/shallow-water/qfunctions/advection.h new file mode 100644 index 0000000000..59bdb26480 --- /dev/null +++ b/examples/fluids/shallow-water/qfunctions/advection.h @@ -0,0 +1,433 @@ +// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at +// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights +// reserved. See files LICENSE and NOTICE for details. +// +// This file is part of CEED, a collection of benchmarks, miniapps, software +// libraries and APIs for efficient high-order finite element and spectral +// element discretizations for exascale applications. For more information and +// source code availability see http://github.com/ceed. +// +// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, +// a collaborative effort of two U.S. Department of Energy organizations (Office +// of Science and the National Nuclear Security Administration) responsible for +// the planning and preparation of a capable exascale ecosystem, including +// software, applications, hardware, advanced system engineering and early +// testbed platforms, in support of the nation's exascale computing imperative. + +/// @file +/// Initial condition and operator for the shallow-water example using PETSc + +#ifndef advection_h +#define advection_h + +#include "../sw_headers.h" + +#ifndef __CUDACC__ +# include +#endif + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +// ***************************************************************************** +// This QFunction sets the the initial condition for the advection of a cosine +// bell shaped scalar function (test case 1 in "A Standard Test Set for +// Numerical Approximations to the Shallow Water Equations in Spherical +// Geometry" by Williamson et al. (1992) +// ***************************************************************************** +static inline int Exact_SW_Advection(CeedInt dim, CeedScalar time, + const CeedScalar X[], CeedInt Nf, + CeedScalar q[], void *ctx) { + + // Context + const PhysicsContext context = (PhysicsContext)ctx; + const CeedScalar u0 = context->u0; + const CeedScalar h0 = context->h0; + const CeedScalar R = context->R; + const CeedScalar gamma = context->gamma; + const CeedScalar rho = R / 3.; + + // Setup + // -- Compute latitude + const CeedScalar theta = asin(X[2] / R); + // -- Compute longitude + const CeedScalar lambda = atan2(X[1], X[0]); + // -- Compute great circle distance between (lambda, theta) and the center, + // (lambda_c, theta_c) + const CeedScalar lambda_c = 3. * M_PI / 2.; + const CeedScalar theta_c = 0.; + const CeedScalar r = R * acos(sin(theta_c)*sin(theta) + cos(theta_c)*cos(theta)*cos(lambda-lambda_c)); + + // Initial Conditions + q[0] = u0 * (cos(theta)*cos(gamma) + sin(theta)*cos(lambda)*sin(gamma)); + q[1] = -u0 * sin(lambda)*sin(gamma); + q[2] = r < rho ? .5*h0*(1 + cos(M_PI*r/rho)) : 0.; // cosine bell + // Return + return 0; +} + +// ***************************************************************************** +// Initial conditions +// ***************************************************************************** +CEED_QFUNCTION(ICsSW_Advection)(void *ctx, CeedInt Q, + const CeedScalar *const *in, CeedScalar *const *out) { + // Inputs + const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; + + // Outputs + CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; + + CeedPragmaSIMD + // Quadrature Point Loop + for (CeedInt i=0; iH0; + const CeedScalar CtauS = context->CtauS; + const CeedScalar strong_form = context->strong_form; + + CeedPragmaSIMD + // Quadrature Point Loop + for (CeedInt i=0; iH0; + const CeedScalar CtauS = context->CtauS; + const CeedScalar strong_form = context->strong_form; + + CeedPragmaSIMD + // Quadrature Point Loop + for (CeedInt i=0; istabilization) { + case 0: + break; + case 1: dv[j][2][i] += wdetJ * TauS * strongConv * uX[j]; //SU + break; + case 2: dv[j][2][i] += wdetJ * TauS * strongResid * uX[j]; //SUPG + break; + } + } // End Quadrature Point Loop + + return 0; + } + +// ***************************************************************************** +// This QFunction implements the Jacobian of of the advection test for the +// shallow-water equations solver. +// +// For this simple test, the equation represents the tranport of the scalar +// field h advected by the wind u, on a spherical surface, where the state +// variables, u_lambda, u_theta (or u_1, u_2) represent the longitudinal +// and latitudinal components of the velocity (wind) field, and h, represents +// the height function. +// +// Discrete Jacobian: +// dF/dq^n = sigma * dF/dqdot|q^n + dF/dq|q^n +// ("sigma * dF/dqdot|q^n" will be added later) +// ***************************************************************************** +CEED_QFUNCTION(SWJacobian_Advection)(void *ctx, CeedInt Q, + const CeedScalar *const *in, + CeedScalar *const *out) { + // *INDENT-OFF* + // Inputs + const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], + (*deltaq)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[1], + (*qdata)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; + // Outputs + CeedScalar (*deltadvdX)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[1]; + // *INDENT-ON* + // Context + ProblemContext context = (ProblemContext)ctx; + const CeedScalar H0 = context->H0; + + CeedPragmaSIMD + // Quadrature Point Loop + for (CeedInt i=0; i @@ -28,51 +30,44 @@ #define M_PI 3.14159265358979323846 #endif -#ifndef physics_context_struct -#define physics_context_struct -typedef struct { - CeedScalar u0; - CeedScalar v0; - CeedScalar h0; - CeedScalar Omega; - CeedScalar R; - CeedScalar g; - CeedScalar H0; - CeedScalar time; -} PhysicsContext_s; -typedef PhysicsContext_s *PhysicsContext; -#endif // physics_context_struct - // ***************************************************************************** -// This QFunction sets the the initial conditions and boundary conditions -// -// For now we have sinusoidal terrain and constant reference height H0 -// +// This QFunction sets the the initial condition for the steady state nonlinear +// zonal geostrophic flow (test case 2 in "A Standard Test Set for Numerical +// Approximations to the Shallow Water Equations in Spherical Geometry" +// by Williamson et al. (1992) // ***************************************************************************** -static inline int Exact_SW(CeedInt dim, CeedScalar time, const CeedScalar X[], +static inline int Exact_SW(CeedInt dim, CeedScalar time, const CeedScalar X[], CeedInt Nf, CeedScalar q[], void *ctx) { // Context const PhysicsContext context = (PhysicsContext)ctx; const CeedScalar u0 = context->u0; - const CeedScalar v0 = context->v0; const CeedScalar h0 = context->h0; + const CeedScalar R = context->R; + const CeedScalar gamma = context->gamma; + const CeedScalar rho = R / 3.; // Setup - // -- Coordinates - const CeedScalar x = X[0]; - const CeedScalar y = X[1]; + // -- Compute latitude + const CeedScalar theta = asin(X[2] / R); + // -- Compute longitude + const CeedScalar lambda = atan2(X[1], X[0]); + // -- Compute great circle distance between (lambda, theta) and the center, + // (lambda_c, theta_c) + const CeedScalar lambda_c = 3. * M_PI / 2.; + const CeedScalar theta_c = 0.; + const CeedScalar r = R * acos(sin(theta_c)*sin(theta) + cos(theta_c)*cos(theta)*cos(lambda-lambda_c)); // Initial Conditions - q[0] = u0; - q[1] = v0; - q[2] = h0; + q[0] = u0 * (cos(theta)*cos(gamma) + sin(theta)*cos(lambda)*sin(gamma)); + q[1] = -u0 * sin(lambda)*sin(gamma); + q[2] = r < rho ? .5*h0*(1 + cos(M_PI*r/rho)) : 0.; // cosine bell // Return return 0; } // ***************************************************************************** -// Initial conditions for shallow-water +// Initial conditions // ***************************************************************************** CEED_QFUNCTION(ICsSW)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { @@ -98,7 +93,6 @@ CEED_QFUNCTION(ICsSW)(void *ctx, CeedInt Q, return 0; } - // ***************************************************************************** // This QFunction implements the explicit terms of the shallow-water // equations @@ -121,11 +115,11 @@ CEED_QFUNCTION(SWExplicit)(void *ctx, CeedInt Q, const CeedScalar *const *in, // Inputs const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1], - (*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], + (*dq)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[1], (*qdata)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; // Outputs CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], - (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; + (*dv)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[1]; // *INDENT-ON* // Context @@ -209,8 +203,8 @@ CEED_QFUNCTION(SWImplicit)(void *ctx, CeedInt Q, const CeedScalar *const *in, // *INDENT-OFF* // Inputs const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], - (*qdot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1], - (*qdata)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; + (*qdot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], + (*qdata)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; // Outputs CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], (*dv)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[1]; @@ -299,7 +293,7 @@ CEED_QFUNCTION(SWJacobian)(void *ctx, CeedInt Q, const CeedScalar *const *in, // *INDENT-OFF* // Inputs const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], - (*deltaq)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[1], + (*deltaq)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[1], (*qdata)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; // Outputs CeedScalar (*deltadvdX)[3][CEED_Q_VLA] = (CeedScalar(*)[3][CEED_Q_VLA])out[1]; @@ -307,6 +301,7 @@ CEED_QFUNCTION(SWJacobian)(void *ctx, CeedInt Q, const CeedScalar *const *in, // Context const PhysicsContext context = (PhysicsContext)ctx; const CeedScalar g = context->g; + const CeedScalar H0 = context->H0; CeedPragmaSIMD // Quadrature Point Loop @@ -314,8 +309,6 @@ CEED_QFUNCTION(SWJacobian)(void *ctx, CeedInt Q, const CeedScalar *const *in, // Setup // Interp in const CeedScalar h = q[2][i]; - // H0 - const CeedScalar H_0 = q[4][i]; // Functional derivatives in const CeedScalar deltau[2] = {deltaq[0][0][i], deltaq[1][0][i] @@ -347,8 +340,8 @@ CEED_QFUNCTION(SWJacobian)(void *ctx, CeedInt Q, const CeedScalar *const *in, deltadvdX[1][1][i] = - g*wdetJ*dXdx[1][1]*deltah; // theta component // Jacobian spatial terms for F_3(t,q): // - dv \cdot ((H_0 + h) delta u) - deltadvdX[1][2][i] = - (H_0 + h)*wdetJ*(deltau[0]*dXdx[1][0] + deltau[1]*dXdx[1][1]); // theta component - deltadvdX[0][2][i] = - (H_0 + h)*wdetJ*(deltau[0]*dXdx[0][0] + deltau[1]*dXdx[0][1]); // lambda component + deltadvdX[1][2][i] = - (H0 + h)*wdetJ*(deltau[0]*dXdx[1][0] + deltau[1]*dXdx[1][1]); // lambda component + deltadvdX[0][2][i] = - (H0 + h)*wdetJ*(deltau[0]*dXdx[0][0] + deltau[1]*dXdx[0][1]); // theta component } // End Quadrature Point Loop @@ -358,4 +351,4 @@ CEED_QFUNCTION(SWJacobian)(void *ctx, CeedInt Q, const CeedScalar *const *in, // ***************************************************************************** -#endif +#endif // shallowwater_h diff --git a/examples/fluids/shallow-water/shallowwater.c b/examples/fluids/shallow-water/shallowwater.c index 057dc7fb19..4853e725f8 100644 --- a/examples/fluids/shallow-water/shallowwater.c +++ b/examples/fluids/shallow-water/shallowwater.c @@ -54,6 +54,10 @@ int main(int argc, char **argv) { Mat J, interpviz; User user; Units units; + problemType problemChoice; + problemData *problem = NULL; + StabilizationType stab; + PetscBool implicit; PetscInt degree, qextra, outputfreq, steps, contsteps; PetscMPIInt rank; PetscScalar ftime; @@ -68,12 +72,15 @@ int main(int argc, char **argv) { Ceed ceed; CeedData ceeddata; CeedMemType memtyperequested; - PetscScalar meter = 1e-2; // 1 meter in scaled length units - PetscScalar second = 1e-2; // 1 second in scaled time units - PetscScalar Omega = 7.29212e-5; // Earth rotation rate (1/s) - PetscScalar R = 6.37122e6; // Earth radius (m) - PetscScalar g = 9.81; // gravitational acceleration (m/s^2) - PetscScalar H0 = 0; // constant mean height (m) + CeedScalar CtauS = 0.; // dimensionless + CeedScalar strong_form = 0.; // [0,1] + PetscScalar meter = 1e-2; // 1 meter in scaled length units + PetscScalar second = 1e-2; // 1 second in scaled time units + PetscScalar Omega = 7.29212e-5; // Earth rotation rate (1/s) + PetscScalar R = 6.37122e6; // Earth radius (m) + PetscScalar g = 9.81; // gravitational acceleration (m/s^2) + PetscScalar H0 = 0; // constant mean height (m) + PetscScalar gamma = 0; // angle between axis of rotation and polar axis PetscScalar mpersquareds; // Check PETSc CUDA support PetscBool petschavecuda, setmemtyperequest = PETSC_FALSE; @@ -104,6 +111,38 @@ int main(int argc, char **argv) { sizeof(ceedresource), NULL); CHKERRQ(ierr); ierr = PetscOptionsBool("-test", "Run in test mode", NULL, test=PETSC_FALSE, &test, NULL); CHKERRQ(ierr); + problemChoice = SWE_ADVECTION; + ierr = PetscOptionsEnum("-problem", "Problem to solve", NULL, + problemTypes, (PetscEnum)problemChoice, + (PetscEnum *)&problemChoice, NULL); CHKERRQ(ierr); + problem = &problemOptions[problemChoice]; + ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL, + StabilizationTypes, (PetscEnum)(stab = STAB_NONE), + (PetscEnum *)&stab, NULL); CHKERRQ(ierr); + ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation", + NULL, implicit=PETSC_FALSE, &implicit, NULL); + CHKERRQ(ierr); + if (!implicit && stab != STAB_NONE) { + ierr = PetscPrintf(comm, "Warning! Use -stab only with -implicit\n"); + CHKERRQ(ierr); + } + ierr = PetscOptionsScalar("-CtauS", + "Scale coefficient for tau (nondimensional)", + NULL, CtauS, &CtauS, NULL); CHKERRQ(ierr); + if (stab == STAB_NONE && CtauS != 0) { + ierr = PetscPrintf(comm, + "Warning! Use -CtauS only with -stab su or -stab supg\n"); + CHKERRQ(ierr); + } + ierr = PetscOptionsScalar("-strong_form", + "Strong (1) or weak/integrated by parts (0) advection residual", + NULL, strong_form, &strong_form, NULL); + CHKERRQ(ierr); + if (problemChoice == SWE_GEOSTROPHIC && (CtauS != 0 || strong_form != 0)) { + ierr = PetscPrintf(comm, + "Warning! Problem geostrophic does not support -CtauS or -strong_form\n"); + CHKERRQ(ierr); + } ierr = PetscOptionsInt("-viz_refine", "Regular refinement levels for visualization", NULL, viz_refine, &viz_refine, NULL); @@ -130,6 +169,8 @@ int main(int argc, char **argv) { NULL, g, &g, NULL); CHKERRQ(ierr); ierr = PetscOptionsScalar("-H0", "Mean height", NULL, H0, &H0, NULL); CHKERRQ(ierr); + ierr = PetscOptionsScalar("-H0", "Angle between axis of rotation and polar axis", + NULL, gamma, &gamma, NULL); CHKERRQ(ierr); PetscStrncpy(user->outputfolder, ".", 2); ierr = PetscOptionsString("-of", "Output folder", NULL, user->outputfolder, user->outputfolder, @@ -157,7 +198,7 @@ int main(int argc, char **argv) { g *= mpersquareds; // Set up the libCEED context - PhysicsContext_s ctx = { + PhysicsContext_s phys_ctx = { .u0 = 0., .v0 = 0., .h0 = .1, @@ -165,8 +206,16 @@ int main(int argc, char **argv) { .R = R, .g = g, .H0 = H0, + .gamma = gamma, .time = 0. }; + + ProblemContext_s probl_ctx = { + .H0 = H0, + .CtauS = CtauS, + .strong_form = strong_form, + .stabilization = stab + }; // Setup DM if (read_mesh) { @@ -277,6 +326,9 @@ int main(int argc, char **argv) { ierr = PetscPrintf(comm, "\n-- CEED Shallow-water equations solver on the cubed-sphere -- libCEED + PETSc --\n" " rank(s) : %d\n" + " Problem:\n" + " Problem Name : %s\n" + " Stabilization : %s\n" " PETSc:\n" " PETSc Vec Type : %s\n" " libCEED:\n" @@ -292,7 +344,8 @@ int main(int argc, char **argv) { " DoFs per node : %D\n" " Global DoFs : %D\n" " Owned DoFs : %D\n", - comm_size, ceedresource, usedresource, + comm_size, problemTypes[problemChoice], + StabilizationTypes[stab], ceedresource, usedresource, CeedMemTypes[memtypebackend], (setmemtyperequest) ? CeedMemTypes[memtyperequested] : "none", numP, numQ, @@ -307,8 +360,8 @@ int main(int argc, char **argv) { // Setup libCEED's objects ierr = PetscMalloc1(1, &ceeddata); CHKERRQ(ierr); - ierr = SetupLibceed(dm, ceed, degree, topodim, qextra, - ncompx, ncompq, user, ceeddata, &ctx); CHKERRQ(ierr); + ierr = SetupLibceed(dm, ceed, degree, qextra, ncompx, ncompq, user, ceeddata, + problem, &phys_ctx, &probl_ctx); CHKERRQ(ierr); // Set up PETSc context // Set up units structure @@ -343,7 +396,7 @@ int main(int argc, char **argv) { // Apply IC operator and fix multiplicity of initial state vector ierr = ICs_FixMultiplicity(ceeddata->op_ics, ceeddata->xcorners, user->q0ceed, dm, Qloc, Q, ceeddata->Erestrictq, - &ctx, 0.0); CHKERRQ(ierr); + &phys_ctx, 0.0); CHKERRQ(ierr); MPI_Comm_rank(comm, &rank); if (!rank) { diff --git a/examples/fluids/shallow-water/src/setup.c b/examples/fluids/shallow-water/src/setup.c index 585bac2b19..03572295c8 100644 --- a/examples/fluids/shallow-water/src/setup.c +++ b/examples/fluids/shallow-water/src/setup.c @@ -25,13 +25,47 @@ #include #include "../sw_headers.h" // Function prototytes #include "../qfunctions/setup_geo.h" // Geometric factors -#include "../qfunctions/shallowwater.h" // Physics point-wise functions +#include "../qfunctions/advection.h" // Physics point-wise functions +#include "../qfunctions/geostrophic.h" // Physics point-wise functions #if PETSC_VERSION_LT(3,14,0) # define DMPlexGetClosureIndices(a,b,c,d,e,f,g,h,i) DMPlexGetClosureIndices(a,b,c,d,f,g,i) # define DMPlexRestoreClosureIndices(a,b,c,d,e,f,g,h,i) DMPlexRestoreClosureIndices(a,b,c,d,f,g,i) #endif +problemData problemOptions[] = { + [SWE_ADVECTION] = { + .topodim = 2, + .qdatasize = 11, + .setup = SetupGeo, + .setup_loc = SetupGeo_loc, + .ics = ICsSW_Advection, + .ics_loc = ICsSW_Advection_loc, + .apply_explfunction = SWExplicit_Advection, + .apply_explfunction_loc = SWExplicit_Advection_loc, + .apply_implfunction = SWImplicit_Advection, + .apply_implfunction_loc = SWImplicit_Advection_loc, + .apply_jacobian = SWJacobian_Advection, + .apply_jacobian_loc = SWJacobian_Advection_loc, + .non_zero_time = PETSC_TRUE + }, + [SWE_GEOSTROPHIC] = { + .topodim = 2, + .qdatasize = 11, + .setup = SetupGeo, + .setup_loc = SetupGeo_loc, + .ics = ICsSW, + .ics_loc = ICsSW_loc, + .apply_explfunction = SWExplicit, + .apply_explfunction_loc = SWExplicit_loc, + .apply_implfunction = SWImplicit, + .apply_implfunction_loc = SWImplicit_loc, + .apply_jacobian = SWJacobian, + .apply_jacobian_loc = SWJacobian_loc, + .non_zero_time = PETSC_FALSE + } +}; + // ----------------------------------------------------------------------------- // Auxiliary function to create PETSc FE space for a given degree // ----------------------------------------------------------------------------- @@ -263,11 +297,10 @@ PetscErrorCode VectorPlacePetscVec(CeedVector c, Vec p) { // Auxiliary function to set up libCEED objects for a given degree // ----------------------------------------------------------------------------- -PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, - CeedInt topodim, CeedInt qextra, - PetscInt ncompx, PetscInt ncompq, - User user, CeedData data, - PhysicsContext ctx) { +PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, CeedInt qextra, + PetscInt ncompx, PetscInt ncompq, User user, + CeedData data, problemData *problem, + PhysicsContext phys_ctx, ProblemContext probl_ctx) { int ierr; DM dmcoord; Vec Xloc; @@ -278,7 +311,8 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, CeedOperator op_setup, op_ics, op_explicit, op_implicit, op_jacobian; CeedVector xcorners, qdata, q0ceed; - CeedInt P, Q, cStart, cEnd, nelem, qdatasize = 11; + CeedInt P, Q, cStart, cEnd, nelem, qdatasize = problem->qdatasize, + topodim = problem->topodim;; // CEED bases P = degree + 1; @@ -317,20 +351,21 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, user->q0ceed = q0ceed; // Create the Q-Function that builds the quadrature data - CeedQFunctionCreateInterior(ceed, 1, SetupGeo, SetupGeo_loc, &qf_setup); + CeedQFunctionCreateInterior(ceed, 1, problem->setup, problem->setup_loc, + &qf_setup); CeedQFunctionAddInput(qf_setup, "x", ncompx, CEED_EVAL_INTERP); CeedQFunctionAddInput(qf_setup, "dx", ncompx*topodim, CEED_EVAL_GRAD); CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT); CeedQFunctionAddOutput(qf_setup, "qdata", qdatasize, CEED_EVAL_NONE); // Create the Q-Function that sets the ICs of the operator - CeedQFunctionCreateInterior(ceed, 1, ICsSW, ICsSW_loc, &qf_ics); + CeedQFunctionCreateInterior(ceed, 1, problem->ics, problem->ics_loc, &qf_ics); CeedQFunctionAddInput(qf_ics, "x", ncompx, CEED_EVAL_INTERP); CeedQFunctionAddOutput(qf_ics, "q0", ncompq, CEED_EVAL_NONE); // Create the Q-Function that defines the explicit part of the PDE operator - CeedQFunctionCreateInterior(ceed, 1, SWExplicit, SWExplicit_loc, - &qf_explicit); + CeedQFunctionCreateInterior(ceed, 1, problem->apply_explfunction, + problem->apply_explfunction_loc, &qf_explicit); CeedQFunctionAddInput(qf_explicit, "x", ncompx, CEED_EVAL_INTERP); CeedQFunctionAddInput(qf_explicit, "q", ncompq, CEED_EVAL_INTERP); CeedQFunctionAddInput(qf_explicit, "dq", ncompq*topodim, CEED_EVAL_GRAD); @@ -339,8 +374,8 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, CeedQFunctionAddOutput(qf_explicit, "dv", ncompq*topodim, CEED_EVAL_GRAD); // Create the Q-Function that defines the implicit part of the PDE operator - CeedQFunctionCreateInterior(ceed, 1, SWImplicit, SWImplicit_loc, - &qf_implicit); + CeedQFunctionCreateInterior(ceed, 1, problem->apply_implfunction, + problem->apply_implfunction_loc, &qf_implicit); CeedQFunctionAddInput(qf_implicit, "q", ncompq, CEED_EVAL_INTERP); CeedQFunctionAddInput(qf_implicit, "dq", ncompq*topodim, CEED_EVAL_GRAD); CeedQFunctionAddInput(qf_implicit, "qdot", ncompq, CEED_EVAL_INTERP); @@ -350,8 +385,8 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, CeedQFunctionAddOutput(qf_implicit, "dv", ncompq*topodim, CEED_EVAL_GRAD); // Create the Q-Function that defines the action of the Jacobian operator - CeedQFunctionCreateInterior(ceed, 1, SWJacobian, SWJacobian_loc, - &qf_jacobian); + CeedQFunctionCreateInterior(ceed, 1, problem->apply_jacobian, + problem->apply_jacobian_loc, &qf_jacobian); CeedQFunctionAddInput(qf_jacobian, "q", 3, CEED_EVAL_INTERP); CeedQFunctionAddInput(qf_jacobian, "deltaq", 3, CEED_EVAL_NONE); CeedQFunctionAddInput(qf_jacobian, "qdata", 10, CEED_EVAL_NONE); @@ -422,10 +457,10 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, user->op_jacobian = op_jacobian; // Set up the libCEED context - CeedQFunctionSetContext(qf_ics, ctx, sizeof *ctx); - CeedQFunctionSetContext(qf_explicit, ctx, sizeof *ctx); - CeedQFunctionSetContext(qf_implicit, ctx, sizeof *ctx); - CeedQFunctionSetContext(qf_jacobian, ctx, sizeof *ctx); + CeedQFunctionSetContext(qf_ics, phys_ctx, sizeof *phys_ctx); + CeedQFunctionSetContext(qf_explicit, probl_ctx, sizeof *probl_ctx); + CeedQFunctionSetContext(qf_implicit, probl_ctx, sizeof *probl_ctx); + CeedQFunctionSetContext(qf_jacobian, probl_ctx, sizeof *probl_ctx); // Save libCEED data required for level // TODO: check how many of these are really needed outside data->basisx = basisx; diff --git a/examples/fluids/shallow-water/sw_headers.h b/examples/fluids/shallow-water/sw_headers.h index 01e39c93f3..6b216b5c88 100644 --- a/examples/fluids/shallow-water/sw_headers.h +++ b/examples/fluids/shallow-water/sw_headers.h @@ -24,8 +24,10 @@ #include #include -#ifndef physics_context_struct -#define physics_context_struct +// ----------------------------------------------------------------------------- +// Data Structs +// ----------------------------------------------------------------------------- + typedef struct { CeedScalar u0; CeedScalar v0; @@ -35,26 +37,58 @@ typedef struct { CeedScalar g; CeedScalar H0; CeedScalar time; + CeedScalar gamma; } PhysicsContext_s; typedef PhysicsContext_s *PhysicsContext; -#endif // physics_context_struct -// MemType Options -static const char *const memTypes[] = { - "host", - "device", - "memType", "CEED_MEM_", NULL -}; +typedef struct { + CeedScalar H0; + CeedScalar CtauS; + CeedScalar strong_form; + int stabilization; // See StabilizationType: 0=none, 1=SU, 2=SUPG +} ProblemContext_s; +typedef ProblemContext_s *ProblemContext; // Problem specific data typedef struct { CeedInt topodim, qdatasize; - CeedQFunctionUser setup, ics, apply_explfunction, apply_implfunction; + CeedQFunctionUser setup, ics, apply_explfunction, apply_implfunction, + apply_jacobian; const char *setup_loc, *ics_loc, *apply_explfunction_loc, - *apply_implfunction_loc; + *apply_implfunction_loc, *apply_jacobian_loc; const bool non_zero_time; } problemData; +// MemType Options +static const char *const memTypes[] = { + "host", + "device", + "memType", "CEED_MEM_", NULL +}; + +// Problem Options +typedef enum { + SWE_ADVECTION = 0, + SWE_GEOSTROPHIC = 1 +} problemType; +static const char *const problemTypes[] = { + "advection", + "geostrophic", + "problemType", "SWE_", NULL +}; + +typedef enum { + STAB_NONE = 0, + STAB_SU = 1, // Streamline Upwind + STAB_SUPG = 2, // Streamline Upwind Petrov-Galerkin +} StabilizationType; +static const char *const StabilizationTypes[] = { + "none", + "SU", + "SUPG", + "StabilizationType", "STAB_", NULL +}; + // PETSc user data typedef struct User_ *User; typedef struct Units_ *Units; @@ -82,10 +116,6 @@ struct Units_ { PetscScalar mpersquareds; }; -// ----------------------------------------------------------------------------- -// libCEED Data Struct -// ----------------------------------------------------------------------------- - // libCEED data struct typedef struct CeedData_ *CeedData; struct CeedData_ { @@ -98,6 +128,9 @@ struct CeedData_ { CeedVector xcorners, xceed, qdata, q0ceed, mceed, hsceed, H0ceed; }; +// External variables +extern problemData problemOptions[]; + // ----------------------------------------------------------------------------- // Setup DM functions // ----------------------------------------------------------------------------- @@ -119,11 +152,10 @@ PetscErrorCode CreateRestrictionPlex(Ceed ceed, DM dm, CeedInt P, CeedInt ncomp, CeedElemRestriction *Erestrict); // Auxiliary function to set up libCEED objects for a given degree -PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, - CeedInt topodim, CeedInt qextra, - PetscInt ncompx, PetscInt ncompq, - User user, CeedData data, - PhysicsContext ctx); +PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, CeedInt qextra, + PetscInt ncompx, PetscInt ncompq, User user, + CeedData data, problemData *problem, + PhysicsContext phys_ctx, ProblemContext probl_ctx); // ----------------------------------------------------------------------------- // RHS (Explicit part in time-stepper) function setup