Skip to content

Commit

Permalink
Extract generic Navier-Stokes code into 'newtonian.h' (#903)
Browse files Browse the repository at this point in the history 
* feat(fluids): Extract newtownian.h from DC problem

 - Also trade out `int` for `StabilizationType` enum

* style: Fix parameter spacing

* feat: Move generic setup from DC -> newtonian

 - Includes moving from `DCContext` to `NewtonianIdealGasContext` at
   appropriate locations

* feat: Add "still" IC for generic Newtonian IG problems

* feat: Move SetupContext from DC -> newtonian.c

* fix: Remove unused variables, general cleanup

* examples/fluids: remove unnecessary includes

* examples/fluids: copy declarations to make newtonian.h self-contained

* examples/fluids: full name for context

Co-authored-by: Jed Brown <[email protected]>
  • Loading branch information
@jrwrigh @jedbrown
jrwrigh and jedbrown authored Mar 5, 2022
1 parent 4e2ba0c commit 88b783a
Show file tree
Hide file tree
Showing 8 changed files with 1,000 additions and 796 deletions.
4 changes: 2 additions & 2 deletions examples/fluids/navierstokes.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -305,7 +305,7 @@ int main(int argc, char **argv) {

// -- Contexts
CeedQFunctionContextDestroy(&ceed_data->setup_context);
CeedQFunctionContextDestroy(&ceed_data->dc_context);
CeedQFunctionContextDestroy(&ceed_data->newt_ig_context);
CeedQFunctionContextDestroy(&ceed_data->advection_context);
CeedQFunctionContextDestroy(&ceed_data->euler_context);

Expand Down Expand Up @@ -367,7 +367,7 @@ int main(int argc, char **argv) {
ierr = PetscFree(problem); CHKERRQ(ierr);
ierr = PetscFree(bc); CHKERRQ(ierr);
ierr = PetscFree(setup_ctx); CHKERRQ(ierr);
ierr = PetscFree(phys_ctx->dc_ctx); CHKERRQ(ierr);
ierr = PetscFree(phys_ctx->newtonian_ig_ctx); CHKERRQ(ierr);
ierr = PetscFree(phys_ctx->euler_ctx); CHKERRQ(ierr);
ierr = PetscFree(phys_ctx->advection_ctx); CHKERRQ(ierr);
ierr = PetscFree(phys_ctx); CHKERRQ(ierr);
Expand Down
45 changes: 33 additions & 12 deletions examples/fluids/navierstokes.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,7 @@ struct AppCtx_private {
// libCEED data struct
struct CeedData_private {
CeedVector x_coord, q_data;
CeedQFunctionContext setup_context, dc_context, advection_context,
CeedQFunctionContext setup_context, newt_ig_context, advection_context,
euler_context;
CeedQFunction qf_setup_vol, qf_ics, qf_rhs_vol, qf_ifunction_vol,
qf_setup_sur, qf_apply_inflow, qf_apply_outflow;
Expand Down Expand Up @@ -262,19 +262,35 @@ struct AdvectionContext_ {
};
#endif

// Newtonian Ideal Gas
#ifndef newtonian_context_struct
#define newtonian_context_struct
typedef struct NewtonianIdealGasContext_ *NewtonianIdealGasContext;
struct NewtonianIdealGasContext_ {
CeedScalar lambda;
CeedScalar mu;
CeedScalar k;
CeedScalar cv;
CeedScalar cp;
CeedScalar g;
CeedScalar c_tau;
StabilizationType stabilization;
};
#endif

// Struct that contains all enums and structs used for the physics of all problems
struct Physics_private {
DCContext dc_ctx;
EulerContext euler_ctx;
AdvectionContext advection_ctx;
WindType wind_type;
BubbleType bubble_type;
BubbleContinuityType bubble_continuity_type;
EulerTestType euler_test;
StabilizationType stab;
PetscBool implicit;
PetscBool has_curr_time;
PetscBool has_neumann;
NewtonianIdealGasContext newtonian_ig_ctx;
EulerContext euler_ctx;
AdvectionContext advection_ctx;
WindType wind_type;
BubbleType bubble_type;
BubbleContinuityType bubble_continuity_type;
EulerTestType euler_test;
StabilizationType stab;
PetscBool implicit;
PetscBool has_curr_time;
PetscBool has_neumann;
};

// Problem specific data
Expand All @@ -298,6 +314,8 @@ typedef struct {
// Set up problems
// -----------------------------------------------------------------------------
// Set up function for each problem
extern PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm,
void *setup_ctx, void *ctx);
extern PetscErrorCode NS_DENSITY_CURRENT(ProblemData *problem, DM dm,
void *setup_ctx, void *ctx);
extern PetscErrorCode NS_EULER_VORTEX(ProblemData *problem, DM dm,
Expand All @@ -308,6 +326,9 @@ extern PetscErrorCode NS_ADVECTION2D(ProblemData *problem, DM dm,
void *setup_ctx, void *ctx);

// Set up context for each problem
extern PetscErrorCode SetupContext_NEWTONIAN_IG(Ceed ceed, CeedData ceed_data,
AppCtx app_ctx, SetupContext setup_ctx, Physics phys);

extern PetscErrorCode SetupContext_DENSITY_CURRENT(Ceed ceed,
CeedData ceed_data, AppCtx app_ctx, SetupContext setup_ctx, Physics phys);

Expand Down
186 changes: 14 additions & 172 deletions examples/fluids/problems/densitycurrent.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,101 +18,39 @@
/// Utility functions for setting up DENSITY_CURRENT

#include "../navierstokes.h"
#include "../qfunctions/setupgeo.h"
#include "../qfunctions/densitycurrent.h"

PetscErrorCode NS_DENSITY_CURRENT(ProblemData *problem, DM dm, void *setup_ctx,
void *ctx) {
SetupContext setup_context = *(SetupContext *)setup_ctx;
User user = *(User *)ctx;
StabilizationType stab;
MPI_Comm comm = PETSC_COMM_WORLD;
PetscBool implicit;
PetscBool has_curr_time = PETSC_FALSE;
PetscInt ierr;
PetscFunctionBeginUser;

ierr = PetscCalloc1(1, &user->phys->dc_ctx); CHKERRQ(ierr);
PetscInt ierr;
ierr = NS_NEWTONIAN_IG(problem, dm, setup_ctx, ctx); CHKERRQ(ierr);
SetupContext setup_context = *(SetupContext *)setup_ctx;
User user = *(User *)ctx;
MPI_Comm comm = PETSC_COMM_WORLD;
PetscFunctionBeginUser;

// ------------------------------------------------------
// SET UP DENSITY_CURRENT
// ------------------------------------------------------
problem->dim = 3;
problem->q_data_size_vol = 10;
problem->q_data_size_sur = 4;
problem->setup_vol = Setup;
problem->setup_vol_loc = Setup_loc;
problem->setup_sur = SetupBoundary;
problem->setup_sur_loc = SetupBoundary_loc;
problem->ics = ICsDC;
problem->ics_loc = ICsDC_loc;
problem->apply_vol_rhs = DC;
problem->apply_vol_rhs_loc = DC_loc;
problem->apply_vol_ifunction = IFunction_DC;
problem->apply_vol_ifunction_loc = IFunction_DC_loc;
problem->bc = Exact_DC;
problem->setup_ctx = SetupContext_DENSITY_CURRENT;
problem->non_zero_time = PETSC_FALSE;
problem->print_info = PRINT_DENSITY_CURRENT;
problem->ics = ICsDC;
problem->ics_loc = ICsDC_loc;
problem->bc = Exact_DC;

// ------------------------------------------------------
// Create the libCEED context
// ------------------------------------------------------
CeedScalar theta0 = 300.; // K
CeedScalar thetaC = -15.; // K
CeedScalar P0 = 1.e5; // Pa
CeedScalar N = 0.01; // 1/s
CeedScalar cv = 717.; // J/(kg K)
CeedScalar cp = 1004.; // J/(kg K)
CeedScalar g = 9.81; // m/s^2
CeedScalar lambda = -2./3.; // -
CeedScalar mu = 75.; // Pa s, dynamic viscosity
// mu = 75 is not physical for air, but is good for numerical stability
CeedScalar k = 0.02638; // W/(m K)
CeedScalar c_tau = 0.5; // -
// c_tau = 0.5 is reported as "optimal" in Hughes et al 2010
CeedScalar rc = 1000.; // m (Radius of bubble)
PetscReal center[3], dc_axis[3] = {0, 0, 0};
PetscReal domain_min[3], domain_max[3], domain_size[3];
ierr = DMGetBoundingBox(dm, domain_min, domain_max); CHKERRQ(ierr);
for (int i=0; i<3; i++) domain_size[i] = domain_max[i] - domain_min[i];

// ------------------------------------------------------
// Create the PETSc context
// ------------------------------------------------------
PetscScalar meter = 1e-2; // 1 meter in scaled length units
PetscScalar kilogram = 1e-6; // 1 kilogram in scaled mass units
PetscScalar second = 1e-2; // 1 second in scaled time units
PetscScalar Kelvin = 1; // 1 Kelvin in scaled temperature units
PetscScalar W_per_m_K, Pascal, J_per_kg_K, m_per_squared_s;

// ------------------------------------------------------
// Command line Options
// ------------------------------------------------------
ierr = PetscOptionsBegin(comm, NULL, "Options for DENSITY_CURRENT problem",
NULL); CHKERRQ(ierr);
// -- Physics
ierr = PetscOptionsScalar("-theta0", "Reference potential temperature",
NULL, theta0, &theta0, NULL); CHKERRQ(ierr);
ierr = PetscOptionsScalar("-thetaC", "Perturbation of potential temperature",
NULL, thetaC, &thetaC, NULL); CHKERRQ(ierr);
ierr = PetscOptionsScalar("-P0", "Atmospheric pressure",
NULL, P0, &P0, NULL); CHKERRQ(ierr);
ierr = PetscOptionsScalar("-N", "Brunt-Vaisala frequency",
NULL, N, &N, NULL); CHKERRQ(ierr);
ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume",
NULL, cv, &cv, NULL); CHKERRQ(ierr);
ierr = PetscOptionsScalar("-cp", "Heat capacity at constant pressure",
NULL, cp, &cp, NULL); CHKERRQ(ierr);
ierr = PetscOptionsScalar("-g", "Gravitational acceleration",
NULL, g, &g, NULL); CHKERRQ(ierr);
ierr = PetscOptionsScalar("-lambda",
"Stokes hypothesis second viscosity coefficient",
NULL, lambda, &lambda, NULL); CHKERRQ(ierr);
ierr = PetscOptionsScalar("-mu", "Shear dynamic viscosity coefficient",
NULL, mu, &mu, NULL); CHKERRQ(ierr);
ierr = PetscOptionsScalar("-k", "Thermal conductivity",
NULL, k, &k, NULL); CHKERRQ(ierr);
ierr = PetscOptionsScalar("-rc", "Characteristic radius of thermal bubble",
NULL, rc, &rc, NULL); CHKERRQ(ierr);
for (int i=0; i<3; i++) center[i] = .5*domain_size[i];
Expand All @@ -130,126 +68,30 @@ PetscErrorCode NS_DENSITY_CURRENT(ProblemData *problem, DM dm, void *setup_ctx,
for (int i=0; i<3; i++) dc_axis[i] /= norm;
}
}
ierr = PetscOptionsEnum("-stab", "Stabilization method", NULL,
StabilizationTypes, (PetscEnum)(stab = STAB_NONE),
(PetscEnum *)&stab, NULL); CHKERRQ(ierr);
ierr = PetscOptionsScalar("-c_tau", "Stabilization constant",
NULL, c_tau, &c_tau, NULL); CHKERRQ(ierr);
ierr = PetscOptionsBool("-implicit", "Use implicit (IFunction) formulation",
NULL, implicit=PETSC_FALSE, &implicit, NULL);
CHKERRQ(ierr);

// -- Units
ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units",
NULL, meter, &meter, NULL); CHKERRQ(ierr);
meter = fabs(meter);
ierr = PetscOptionsScalar("-units_kilogram","1 kilogram in scaled mass units",
NULL, kilogram, &kilogram, NULL); CHKERRQ(ierr);
kilogram = fabs(kilogram);
ierr = PetscOptionsScalar("-units_second","1 second in scaled time units",
NULL, second, &second, NULL); CHKERRQ(ierr);
second = fabs(second);
ierr = PetscOptionsScalar("-units_Kelvin",
"1 Kelvin in scaled temperature units",
NULL, Kelvin, &Kelvin, NULL); CHKERRQ(ierr);
Kelvin = fabs(Kelvin);

// -- Warnings
if (stab == STAB_SUPG && !implicit) {
ierr = PetscPrintf(comm,
"Warning! Use -stab supg only with -implicit\n");
CHKERRQ(ierr);
}

ierr = PetscOptionsEnd(); CHKERRQ(ierr);

// ------------------------------------------------------
// Set up the PETSc context
// ------------------------------------------------------
// -- Define derived units
Pascal = kilogram / (meter * PetscSqr(second));
J_per_kg_K = PetscSqr(meter) / (PetscSqr(second) * Kelvin);
m_per_squared_s = meter / PetscSqr(second);
W_per_m_K = kilogram * meter / (pow(second,3) * Kelvin);

user->units->meter = meter;
user->units->kilogram = kilogram;
user->units->second = second;
user->units->Kelvin = Kelvin;
user->units->Pascal = Pascal;
user->units->J_per_kg_K = J_per_kg_K;
user->units->m_per_squared_s = m_per_squared_s;
user->units->W_per_m_K = W_per_m_K;

// ------------------------------------------------------
// Set up the libCEED context
// ------------------------------------------------------
// -- Scale variables to desired units
theta0 *= Kelvin;
thetaC *= Kelvin;
P0 *= Pascal;
N *= (1./second);
cv *= J_per_kg_K;
cp *= J_per_kg_K;
g *= m_per_squared_s;
mu *= Pascal * second;
k *= W_per_m_K;
rc = fabs(rc) * meter;
for (int i=0; i<3; i++) domain_size[i] *= meter;
PetscScalar meter = user->units->meter;
rc = fabs(rc) * meter;
for (int i=0; i<3; i++) center[i] *= meter;
problem->dm_scale = meter;

// -- Setup Context
setup_context->theta0 = theta0;
setup_context->thetaC = thetaC;
setup_context->P0 = P0;
setup_context->N = N;
setup_context->cv = cv;
setup_context->cp = cp;
setup_context->g = g;
setup_context->rc = rc;
setup_context->lx = domain_size[0];
setup_context->ly = domain_size[1];
setup_context->lz = domain_size[2];
setup_context->center[0] = center[0];
setup_context->center[1] = center[1];
setup_context->center[2] = center[2];
setup_context->dc_axis[0] = dc_axis[0];
setup_context->dc_axis[1] = dc_axis[1];
setup_context->dc_axis[2] = dc_axis[2];
setup_context->time = 0;

// -- QFunction Context
user->phys->stab = stab;
user->phys->implicit = implicit;
user->phys->has_curr_time = has_curr_time;
user->phys->dc_ctx->lambda = lambda;
user->phys->dc_ctx->mu = mu;
user->phys->dc_ctx->k = k;
user->phys->dc_ctx->cv = cv;
user->phys->dc_ctx->cp = cp;
user->phys->dc_ctx->g = g;
user->phys->dc_ctx->c_tau = c_tau;
user->phys->dc_ctx->stabilization = stab;

PetscFunctionReturn(0);
}

PetscErrorCode SetupContext_DENSITY_CURRENT(Ceed ceed, CeedData ceed_data,
AppCtx app_ctx, SetupContext setup_ctx, Physics phys) {
PetscFunctionBeginUser;
CeedQFunctionContextCreate(ceed, &ceed_data->setup_context);
CeedQFunctionContextSetData(ceed_data->setup_context, CEED_MEM_HOST,
CEED_USE_POINTER, sizeof(*setup_ctx), setup_ctx);
CeedQFunctionSetContext(ceed_data->qf_ics, ceed_data->setup_context);
CeedQFunctionContextCreate(ceed, &ceed_data->dc_context);
CeedQFunctionContextSetData(ceed_data->dc_context, CEED_MEM_HOST,
CEED_USE_POINTER,
sizeof(*phys->dc_ctx), phys->dc_ctx);
if (ceed_data->qf_rhs_vol)
CeedQFunctionSetContext(ceed_data->qf_rhs_vol, ceed_data->dc_context);
if (ceed_data->qf_ifunction_vol)
CeedQFunctionSetContext(ceed_data->qf_ifunction_vol, ceed_data->dc_context);
PetscInt ierr = SetupContext_NEWTONIAN_IG(ceed, ceed_data, app_ctx, setup_ctx,
phys);
CHKERRQ(ierr);
PetscFunctionReturn(0);
}

Expand Down
Loading

0 comments on commit 88b783a

Please sign in to comment.