Skip to content

Commit

Permalink
Fluids - Newtonian with Primitive variables (#1011)
Browse files Browse the repository at this point in the history 
* Fluids - Initial commit for Newtonian primitive variables

* Fluids - include A0 (dU/dY) in the output

* Fluids - added ICs for IG in primitive variables

* Fluids - added Jacobean QFunction for primitive variables

* Fluids - added RHS QFunction for primitive variables

* Fluids - fixed compilation errors and warnings

* Fluids - added ICs for density_current with primitive variables

* Fluids - In/OutFlow BCs for channel in primitive variables

* Adding the missing parts after rebasing onto main

* Fluids - Use correct component names for primitive variables

* Fluids - Primitive variables supported only with implicit scheme

* Fluids - drop in/outflow for channel flow and call Exact_Channel_Prim() in Exact_Channel()

* Fluids - Set solver's QFunction data in an if-else statement

* Fluids - style

* Fluids - add a comment to explain why the the gravity body force is excluded from the potential energy.

* Fluids - Exact_Channel return State

* Fluids - density_current: some style and cleanup

* Fluids - DC: refactor & cleanup

* Fluids - Singel QFunction for prim&cons

* Fluids - Use absolute temperature

* Fluids - DC: Fix pressure
  • Loading branch information
@LeilaGhaffari
LeilaGhaffari authored Jul 18, 2022
1 parent 4add99e commit dc805cc
Show file tree
Hide file tree
Showing 10 changed files with 670 additions and 153 deletions.
16 changes: 6 additions & 10 deletions examples/fluids/navierstokes.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -194,6 +194,7 @@ struct Physics_private {
EulerTestType euler_test;
StabilizationType stab;
PetscBool implicit;
PetscBool primitive;
PetscBool has_curr_time;
PetscBool has_neumann;
CeedContextFieldLabel solution_time_label;
Expand Down Expand Up @@ -251,20 +252,15 @@ extern PetscErrorCode NS_ADVECTION2D(ProblemData *problem, DM dm,
void *ctx);

// Print function for each problem
extern PetscErrorCode PRINT_DENSITY_CURRENT(ProblemData *problem,
AppCtx app_ctx);
extern PetscErrorCode PRINT_NEWTONIAN(ProblemData *problem, AppCtx app_ctx);

extern PetscErrorCode PRINT_EULER_VORTEX(ProblemData *problem,
AppCtx app_ctx);
extern PetscErrorCode PRINT_EULER_VORTEX(ProblemData *problem, AppCtx app_ctx);

extern PetscErrorCode PRINT_SHOCKTUBE(ProblemData *problem,
AppCtx app_ctx);
extern PetscErrorCode PRINT_SHOCKTUBE(ProblemData *problem, AppCtx app_ctx);

extern PetscErrorCode PRINT_ADVECTION(ProblemData *problem,
AppCtx app_ctx);
extern PetscErrorCode PRINT_ADVECTION(ProblemData *problem, AppCtx app_ctx);

extern PetscErrorCode PRINT_ADVECTION2D(ProblemData *problem,
AppCtx app_ctx);
extern PetscErrorCode PRINT_ADVECTION2D(ProblemData *problem, AppCtx app_ctx);

// -----------------------------------------------------------------------------
// libCEED functions
Expand Down
14 changes: 8 additions & 6 deletions examples/fluids/problems/channel.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,12 +28,14 @@ PetscErrorCode NS_CHANNEL(ProblemData *problem, DM dm, void *ctx) {
// SET UP Channel
// ------------------------------------------------------
CeedQFunctionContextDestroy(&problem->ics.qfunction_context);
problem->ics.qfunction = ICsChannel;
problem->ics.qfunction_loc = ICsChannel_loc;
problem->apply_inflow.qfunction = Channel_Inflow;
problem->apply_inflow.qfunction_loc = Channel_Inflow_loc;
problem->apply_outflow.qfunction = Channel_Outflow;
problem->apply_outflow.qfunction_loc = Channel_Outflow_loc;
problem->ics.qfunction = ICsChannel;
problem->ics.qfunction_loc = ICsChannel_loc;
if (!user->phys->primitive) {
problem->apply_inflow.qfunction = Channel_Inflow;
problem->apply_inflow.qfunction_loc = Channel_Inflow_loc;
problem->apply_outflow.qfunction = Channel_Outflow;
problem->apply_outflow.qfunction_loc = Channel_Outflow_loc;
}

// -- Command Line Options
CeedScalar umax = 10.; // m/s
Expand Down
61 changes: 35 additions & 26 deletions examples/fluids/problems/densitycurrent.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,21 +14,22 @@

PetscErrorCode NS_DENSITY_CURRENT(ProblemData *problem, DM dm, void *ctx) {

PetscInt ierr;
SetupContext setup_context;
User user = *(User *)ctx;
MPI_Comm comm = PETSC_COMM_WORLD;
PetscInt ierr;
MPI_Comm comm = PETSC_COMM_WORLD;
User user = *(User *)ctx;
DensityCurrentContext dc_ctx;
CeedQFunctionContext density_current_context;
NewtonianIdealGasContext newtonian_ig_ctx;

PetscFunctionBeginUser;
ierr = NS_NEWTONIAN_IG(problem, dm, ctx); CHKERRQ(ierr);
ierr = PetscCalloc1(1, &dc_ctx); CHKERRQ(ierr);
// ------------------------------------------------------
// SET UP DENSITY_CURRENT
// ------------------------------------------------------
problem->ics.qfunction = ICsDC;
CeedQFunctionContextDestroy(&problem->ics.qfunction_context);
problem->ics.qfunction = ICsDC;
problem->ics.qfunction_loc = ICsDC_loc;
problem->bc = Exact_DC;
CeedQFunctionContextGetData(problem->ics.qfunction_context, CEED_MEM_HOST,
&setup_context);

// ------------------------------------------------------
// Create the libCEED context
Expand Down Expand Up @@ -85,10 +86,10 @@ PetscErrorCode NS_DENSITY_CURRENT(ProblemData *problem, DM dm, void *ctx) {

PetscOptionsEnd();

PetscScalar meter = user->units->meter;
PetscScalar second = user->units->second;
PetscScalar Kelvin = user->units->Kelvin;
PetscScalar Pascal = user->units->Pascal;
PetscScalar meter = user->units->meter;
PetscScalar second = user->units->second;
PetscScalar Kelvin = user->units->Kelvin;
PetscScalar Pascal = user->units->Pascal;
rc = fabs(rc) * meter;
theta0 *= Kelvin;
thetaC *= Kelvin;
Expand All @@ -97,21 +98,29 @@ PetscErrorCode NS_DENSITY_CURRENT(ProblemData *problem, DM dm, void *ctx) {
for (PetscInt i = 0; i < 3; i++)
center[i] *= meter;

setup_context->theta0 = theta0;
setup_context->thetaC = thetaC;
setup_context->P0 = P0;
setup_context->N = N;
setup_context->rc = rc;
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];
dc_ctx->theta0 = theta0;
dc_ctx->thetaC = thetaC;
dc_ctx->P0 = P0;
dc_ctx->N = N;
dc_ctx->rc = rc;
dc_ctx->center[0] = center[0];
dc_ctx->center[1] = center[1];
dc_ctx->center[2] = center[2];
dc_ctx->dc_axis[0] = dc_axis[0];
dc_ctx->dc_axis[1] = dc_axis[1];
dc_ctx->dc_axis[2] = dc_axis[2];

problem->bc_ctx =
setup_context; // This is bad, context data should only be accessed via Get/Restore
CeedQFunctionContextRestoreData(problem->ics.qfunction_context, &setup_context);
CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context,
CEED_MEM_HOST, &newtonian_ig_ctx);
dc_ctx->newtonian_ctx = *newtonian_ig_ctx;
CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context,
&newtonian_ig_ctx);
CeedQFunctionContextCreate(user->ceed, &density_current_context);
CeedQFunctionContextSetData(density_current_context, CEED_MEM_HOST,
CEED_USE_POINTER, sizeof(*dc_ctx), dc_ctx);
CeedQFunctionContextSetDataDestroy(density_current_context, CEED_MEM_HOST,
FreeContextPetsc);
problem->ics.qfunction_context = density_current_context;

PetscFunctionReturn(0);
}
81 changes: 51 additions & 30 deletions examples/fluids/problems/newtonian.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,12 +38,14 @@ static PetscErrorCode UnitTests_Newtonian(User user,
NewtonianIdealGasContext gas) {

Units units = user->units;
const CeedScalar eps = 1e-6;
const CeedScalar kg = units->kilogram, m = units->meter, sec = units->second,
const CeedScalar eps = 1e-6;
const CeedScalar kg = units->kilogram,
m = units->meter,
sec = units->second,
Pascal = units->Pascal;

PetscFunctionBeginUser;
const CeedScalar rho = 1.2 * kg / (m*m*m), u = 40 * m/sec;
const CeedScalar rho = 1.2 * kg / (m*m*m),
u = 40 * m/sec;
CeedScalar U[5] = {rho, rho*u, rho *u*1.1, rho *u*1.2, 250e3*Pascal + .5*rho *u*u};
const CeedScalar x[3] = {.1, .2, .3};
State s = StateFromU(gas, U, x);
Expand Down Expand Up @@ -74,7 +76,8 @@ PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm, void *ctx) {
StabilizationType stab;
MPI_Comm comm = PETSC_COMM_WORLD;
PetscBool implicit;
PetscBool has_curr_time = PETSC_FALSE, unit_tests;
PetscBool has_curr_time = PETSC_FALSE,
prim_var, unit_tests;
PetscInt ierr;
NewtonianIdealGasContext newtonian_ig_ctx;
CeedQFunctionContext newtonian_ig_context;
Expand All @@ -92,28 +95,12 @@ PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm, void *ctx) {
problem->jac_data_size_sur = 11;
problem->setup_vol.qfunction = Setup;
problem->setup_vol.qfunction_loc = Setup_loc;
problem->ics.qfunction = ICsNewtonianIG;
problem->ics.qfunction_loc = ICsNewtonianIG_loc;
problem->setup_sur.qfunction = SetupBoundary;
problem->setup_sur.qfunction_loc = SetupBoundary_loc;
problem->apply_vol_rhs.qfunction = RHSFunction_Newtonian;
problem->apply_vol_rhs.qfunction_loc = RHSFunction_Newtonian_loc;
problem->apply_vol_ifunction.qfunction = IFunction_Newtonian;
problem->apply_vol_ifunction.qfunction_loc = IFunction_Newtonian_loc;
problem->apply_vol_ijacobian.qfunction = IJacobian_Newtonian;
problem->apply_vol_ijacobian.qfunction_loc = IJacobian_Newtonian_loc;
problem->apply_inflow.qfunction = BoundaryIntegral;
problem->apply_inflow.qfunction_loc = BoundaryIntegral_loc;
problem->apply_inflow_jacobian.qfunction = BoundaryIntegral_Jacobian;
problem->apply_inflow_jacobian.qfunction_loc = BoundaryIntegral_Jacobian_loc;
problem->apply_outflow.qfunction = PressureOutflow;
problem->apply_outflow.qfunction_loc = PressureOutflow_loc;
problem->apply_outflow_jacobian.qfunction = PressureOutflow_Jacobian;
problem->apply_outflow_jacobian.qfunction_loc = PressureOutflow_Jacobian_loc;
problem->bc = NULL;
problem->bc_ctx = setup_context;
problem->non_zero_time = PETSC_FALSE;
problem->print_info = PRINT_DENSITY_CURRENT;
problem->print_info = PRINT_NEWTONIAN;

// ------------------------------------------------------
// Create the libCEED context
Expand All @@ -125,11 +112,11 @@ PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm, void *ctx) {
CeedScalar mu = 1.8e-5; // Pa s, dynamic viscosity
CeedScalar k = 0.02638; // W/(m K)
CeedScalar c_tau = 0.5; // -
CeedScalar Ctau_t = 1.0; // -
CeedScalar Ctau_v = 36.0; // TODO make function of degree
CeedScalar Ctau_C = 1.0; // TODO make function of degree
CeedScalar Ctau_M = 1.0; // TODO make function of degree
CeedScalar Ctau_E = 1.0; // TODO make function of degree
CeedScalar Ctau_t = 1.0; // -
CeedScalar Ctau_v = 36.0; // TODO make function of degree
CeedScalar Ctau_C = 1.0; // TODO make function of degree
CeedScalar Ctau_M = 1.0; // TODO make function of degree
CeedScalar Ctau_E = 1.0; // TODO make function of degree
PetscReal domain_min[3], domain_max[3], domain_size[3];
ierr = DMGetBoundingBox(dm, domain_min, domain_max); CHKERRQ(ierr);
for (PetscInt i=0; i<3; i++) domain_size[i] = domain_max[i] - domain_min[i];
Expand All @@ -140,14 +127,42 @@ PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm, void *ctx) {
PetscScalar meter = 1; // 1 meter in scaled length units
PetscScalar kilogram = 1; // 1 kilogram in scaled mass units
PetscScalar second = 1; // 1 second in scaled time units
PetscScalar Kelvin = 1; // 1 Kelvin in scaled temperature 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
// ------------------------------------------------------
PetscOptionsBegin(comm, NULL, "Options for Newtonian Ideal Gas based problem",
NULL);
// -- Conservative vs Primitive variables
ierr = PetscOptionsBool("-primitive", "Use primitive variables",
NULL, prim_var=PETSC_FALSE, &prim_var, NULL); CHKERRQ(ierr);
if (prim_var) {
problem->ics.qfunction = ICsNewtonianIG_Prim;
problem->ics.qfunction_loc = ICsNewtonianIG_Prim_loc;
problem->apply_vol_ifunction.qfunction = IFunction_Newtonian_Prim;
problem->apply_vol_ifunction.qfunction_loc = IFunction_Newtonian_Prim_loc;
problem->apply_vol_ijacobian.qfunction = IJacobian_Newtonian_Prim;
problem->apply_vol_ijacobian.qfunction_loc = IJacobian_Newtonian_Prim_loc;
} else {
problem->ics.qfunction = ICsNewtonianIG;
problem->ics.qfunction_loc = ICsNewtonianIG_loc;
problem->apply_vol_rhs.qfunction = RHSFunction_Newtonian;
problem->apply_vol_rhs.qfunction_loc = RHSFunction_Newtonian_loc;
problem->apply_vol_ifunction.qfunction = IFunction_Newtonian;
problem->apply_vol_ifunction.qfunction_loc = IFunction_Newtonian_loc;
problem->apply_vol_ijacobian.qfunction = IJacobian_Newtonian;
problem->apply_vol_ijacobian.qfunction_loc = IJacobian_Newtonian_loc;
problem->apply_inflow.qfunction = BoundaryIntegral;
problem->apply_inflow.qfunction_loc = BoundaryIntegral_loc;
problem->apply_inflow_jacobian.qfunction = BoundaryIntegral_Jacobian;
problem->apply_inflow_jacobian.qfunction_loc = BoundaryIntegral_Jacobian_loc;
problem->apply_outflow.qfunction = PressureOutflow;
problem->apply_outflow.qfunction_loc = PressureOutflow_loc;
problem->apply_outflow_jacobian.qfunction = PressureOutflow_Jacobian;
problem->apply_outflow_jacobian.qfunction_loc = PressureOutflow_Jacobian_loc;
}

// -- Physics
ierr = PetscOptionsScalar("-cv", "Heat capacity at constant volume",
Expand Down Expand Up @@ -208,6 +223,10 @@ PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm, void *ctx) {
"Warning! Use -stab supg only with -implicit\n");
CHKERRQ(ierr);
}
if (prim_var && !implicit) {
SETERRQ(comm, PETSC_ERR_ARG_NULL,
"RHSFunction is not provided for primitive variables (use -primitive only with -implicit)\n");
}
PetscOptionsEnd();

// ------------------------------------------------------
Expand Down Expand Up @@ -252,6 +271,7 @@ PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm, void *ctx) {
// -- Solver Settings
user->phys->stab = stab;
user->phys->implicit = implicit;
user->phys->primitive = prim_var;
user->phys->has_curr_time = has_curr_time;

// -- QFunction Context
Expand All @@ -268,6 +288,7 @@ PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm, void *ctx) {
newtonian_ig_ctx->Ctau_E = Ctau_E;
newtonian_ig_ctx->stabilization = stab;
newtonian_ig_ctx->is_implicit = implicit;
newtonian_ig_ctx->primitive = prim_var;
ierr = PetscArraycpy(newtonian_ig_ctx->g, g, 3); CHKERRQ(ierr);

CeedQFunctionContextCreate(user->ceed, &problem->ics.qfunction_context);
Expand Down Expand Up @@ -311,8 +332,8 @@ PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm, void *ctx) {
PetscFunctionReturn(0);
}

PetscErrorCode PRINT_DENSITY_CURRENT(ProblemData *problem,
AppCtx app_ctx) {
PetscErrorCode PRINT_NEWTONIAN(ProblemData *problem, AppCtx app_ctx) {

MPI_Comm comm = PETSC_COMM_WORLD;
PetscErrorCode ierr;
NewtonianIdealGasContext newtonian_ctx;
Expand Down
Loading

0 comments on commit dc805cc

Please sign in to comment.