Skip to content

Commit

Permalink
Merge pull request CEED#1539 from CEED/jrwrigh/explicit_ns
Browse files Browse the repository at this point in the history 
fluids: Make Newtonian BCs compatible with explicit timestepping, add test
  • Loading branch information
@jrwrigh
jrwrigh authored Apr 2, 2024
2 parents d8d1641 + 4c0e823 commit 780fa44
Show file tree
Hide file tree
Showing 9 changed files with 82 additions and 121 deletions.
1 change: 0 additions & 1 deletion examples/fluids/gaussianwave.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,6 @@ ts:
max_steps: 100
dt: 2e-3
type: alpha
alpha_radius: 0.5
#monitor_solution: cgns:nwave.cgns
#monitor_solution_interval: 10

Expand Down
5 changes: 3 additions & 2 deletions examples/fluids/navierstokes.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,17 +18,18 @@
// ./navierstokes -ceed /cpu/self -options_file gaussianwave.yml
// ./navierstokes -ceed /gpu/cuda -problem advection -degree 1
//
//TESTARGS(name="Gaussian Wave, explicit") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/gaussianwave.yaml -compare_final_state_atol 1e-8 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-gaussianwave-explicit.bin -dm_plex_box_faces 2,2,1 -ts_max_steps 5 -degree 3 -implicit false -ts_type rk -stab su -state_var conservative
//TESTARGS(name="Advection 2D, rotation, explicit, su, consistent mass") -ceed {ceed_resource} -test_type solver -problem advection -degree 3 -dm_plex_box_faces 2,2 -dm_plex_box_lower 0,0 -dm_plex_box_upper 125,125 -bc_wall 1,2,3,4 -wall_comps 4 -units_kilogram 1e-9 -rc 100. -ts_dt 1e-3 -ts_max_steps 10 -stab su -Ctaus 0.5 -mass_ksp_type cg -mass_pc_jacobi_type diagonal -compare_final_state_atol 1e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-adv2d-rotation-explicit-stab-su-consistent-mass.bin
//TESTARGS(name="Advection, skew") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/advection.yaml -ts_max_steps 5 -wind_type translation -wind_translation -0.5547002,0.83205029,0 -advection_ic_type skew -dm_plex_box_faces 2,1,1 -degree 2 -stab supg -stab_tau advdiff_shakib -Ctau_a 4 -ksp_type gmres -compare_final_state_atol 5e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-adv-skew.bin
//TESTARGS(name="Blasius, bc_slip") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/blasius.yaml -ts_max_steps 5 -dm_plex_box_faces 3,20,1 -platemesh_nDelta 10 -platemesh_growth 1.2 -bc_outflow 5 -bc_slip 4 -compare_final_state_atol 2E-11 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius-bc_slip.bin
//TESTARGS(name="Blasius, SGS DataDriven Sequential") -ceed {ceed_resource} -options_file examples/fluids/tests-output/blasius_stgtest.yaml -sgs_model_type data_driven -sgs_model_dd_leakyrelu_alpha 0.3 -sgs_model_dd_parameter_dir examples/fluids/dd_sgs_data -ts_dt 2e-9 -state_var primitive -ksp_rtol 1e-12 -snes_rtol 1e-12 -stg_mean_only -stg_fluctuating_IC -test_type solver -compare_final_state_atol 1e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius-sgs-data-driven.bin -sgs_model_dd_use_fused false
//TESTARGS(name="Advection, rotation, cosine") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/advection.yaml -ts_max_steps 0 -advection_ic_type cosine_hill -dm_plex_box_faces 2,1,1 -compare_final_state_atol 1e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-adv-rotation-cosine.bin
//TESTARGS(name="Gaussian Wave, using MatShell") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/gaussianwave.yaml -compare_final_state_atol 1e-8 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-gaussianwave-shell.bin -dm_plex_box_faces 2,2,1 -ts_max_steps 5 -degree 3 -amat_type shell -pc_type vpbjacobi
//TESTARGS(name="Gaussian Wave, using MatShell") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/gaussianwave.yaml -compare_final_state_atol 1e-8 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-gaussianwave-shell.bin -dm_plex_box_faces 2,2,1 -ts_max_steps 5 -degree 3 -amat_type shell -pc_type vpbjacobi -ts_alpha_radius 0.5
//TESTARGS(name="Taylor-Green Vortex IC") -ceed {ceed_resource} -problem taylor_green -test_type solver -dm_plex_dim 3 -dm_plex_box_faces 6,6,6 -ts_max_steps 0 -compare_final_state_atol 1e-12 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-taylor-green-IC.bin
//TESTARGS(name="Blasius, SGS DataDriven Fused") -ceed {ceed_resource} -options_file examples/fluids/tests-output/blasius_stgtest.yaml -sgs_model_type data_driven -sgs_model_dd_leakyrelu_alpha 0.3 -sgs_model_dd_parameter_dir examples/fluids/dd_sgs_data -ts_dt 2e-9 -state_var primitive -ksp_rtol 1e-12 -snes_rtol 1e-12 -stg_mean_only -stg_fluctuating_IC -test_type solver -compare_final_state_atol 1e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius-sgs-data-driven.bin
//TESTARGS(name="Blasius, Anisotropic Differential Filter") -ceed {ceed_resource} -test_type diff_filter -options_file examples/fluids/tests-output/blasius_test.yaml -compare_final_state_atol 5e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius_diff_filter_aniso_vandriest.bin -diff_filter_monitor -ts_max_steps 0 -state_var primitive -diff_filter_friction_length 1e-5 -diff_filter_wall_damping_function van_driest -diff_filter_ksp_rtol 1e-8 -diff_filter_grid_based_width -diff_filter_width_scaling 1,0.7,1
//TESTARGS(name="Blasius, Isotropic Differential Filter") -ceed {ceed_resource} -test_type diff_filter -options_file examples/fluids/tests-output/blasius_test.yaml -compare_final_state_atol 2e-12 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius_diff_filter_iso.bin -diff_filter_monitor -ts_max_steps 0 -diff_filter_width_scaling 4.2e-5,4.2e-5,4.2e-5 -diff_filter_ksp_atol 1e-14 -diff_filter_ksp_rtol 1e-16
//TESTARGS(name="Gaussian Wave, with IDL") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/gaussianwave.yaml -compare_final_state_atol 2e-11 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-gaussianwave-IDL.bin -dm_plex_box_faces 5,5,1 -ts_max_steps 5 -idl_decay_time 2e-3 -idl_length 0.25 -idl_start 0
//TESTARGS(name="Gaussian Wave, with IDL") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/gaussianwave.yaml -compare_final_state_atol 2e-11 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-gaussianwave-IDL.bin -dm_plex_box_faces 5,5,1 -ts_max_steps 5 -idl_decay_time 2e-3 -idl_length 0.25 -idl_start 0 -ts_alpha_radius 0.5
//TESTARGS(name="Spanwise Turbulence Statistics") -ceed {ceed_resource} -test_type turb_spanstats -options_file examples/fluids/tests-output/stats_test.yaml -compare_final_state_atol 1E-11 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-turb-spanstats-stats.bin
//TESTARGS(name="Blasius") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/tests-output/blasius_test.yaml -compare_final_state_atol 2E-11 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius.bin
//TESTARGS(name="Blasius, STG Inflow") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/tests-output/blasius_stgtest.yaml -compare_final_state_atol 2E-11 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius_STG.bin
Expand Down
54 changes: 25 additions & 29 deletions examples/fluids/qfunctions/bc_freestream.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,20 +17,18 @@
// *****************************************************************************
CEED_QFUNCTION_HELPER int Freestream(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var,
RiemannFluxType flux_type) {
const FreestreamContext context = (FreestreamContext)ctx;
const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
const CeedScalar(*q_data_sur) = in[2];
CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
CeedScalar(*jac_data_sur) = context->newtonian_ctx.is_implicit ? out[1] : NULL;

CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
CeedScalar(*jac_data_sur) = out[1];

const FreestreamContext context = (FreestreamContext)ctx;
const NewtonianIdealGasContext newt_ctx = &context->newtonian_ctx;
const bool is_implicit = newt_ctx->is_implicit;
const CeedScalar zeros[6] = {0.};

CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
State s = StateFromQ(newt_ctx, qi, state_var);
const State s = StateFromQ(newt_ctx, qi, state_var);

CeedScalar wdetJb, norm[3];
QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, norm);
Expand All @@ -49,8 +47,11 @@ CEED_QFUNCTION_HELPER int Freestream(void *ctx, CeedInt Q, const CeedScalar *con
UnpackState_U(flux, Flux);
for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j];

StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur);
StoredValuesPack(Q, i, 5, 6, zeros, jac_data_sur); // Every output value must be set
if (is_implicit) {
CeedScalar zeros[6] = {0.};
StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur);
StoredValuesPack(Q, i, 5, 6, zeros, jac_data_sur); // Every output value must be set
}
}
return 0;
}
Expand Down Expand Up @@ -154,16 +155,13 @@ CEED_QFUNCTION_HELPER CeedScalar Softplus_fwd(CeedScalar x, CeedScalar dx, CeedS
// little or no benefit in the tests we've run thus far, thus we recommend
// skipping this feature and just allowing recirculation.
CEED_QFUNCTION_HELPER int RiemannOutflow(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
// Inputs
const OutflowContext outflow = (OutflowContext)ctx;
const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
const CeedScalar(*Grad_q) = in[1];
const CeedScalar(*q_data_sur) = in[2];
CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
CeedScalar(*jac_data_sur) = outflow->gas.is_implicit ? out[1] : NULL;

// Outputs
CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
CeedScalar(*jac_data_sur) = out[1];

const OutflowContext outflow = (OutflowContext)ctx;
const NewtonianIdealGasContext gas = &outflow->gas;
const bool is_implicit = gas->is_implicit;

Expand All @@ -172,7 +170,7 @@ CEED_QFUNCTION_HELPER int RiemannOutflow(void *ctx, CeedInt Q, const CeedScalar
QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, norm);
wdetJb *= is_implicit ? -1. : 1.;
const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
State s_int = StateFromQ(gas, qi, state_var);
const State s_int = StateFromQ(gas, qi, state_var);

StatePrimitive y_ext = s_int.Y;
y_ext.pressure = outflow->pressure;
Expand Down Expand Up @@ -201,8 +199,10 @@ CEED_QFUNCTION_HELPER int RiemannOutflow(void *ctx, CeedInt Q, const CeedScalar
for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j];

// Save values for Jacobian
StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur);
StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur);
if (is_implicit) {
StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur);
StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur);
}
}
return 0;
}
Expand Down Expand Up @@ -299,23 +299,17 @@ CEED_QFUNCTION(RiemannOutflow_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedSca
// acoustics, vortices, etc.
// *****************************************************************************
CEED_QFUNCTION_HELPER int PressureOutflow(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
// Inputs
const OutflowContext outflow = (OutflowContext)ctx;
const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
const CeedScalar(*Grad_q) = in[1];
const CeedScalar(*q_data_sur) = in[2];
CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
CeedScalar(*jac_data_sur) = outflow->gas.is_implicit ? out[1] : NULL;

// Outputs
CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
CeedScalar(*jac_data_sur) = out[1];

const OutflowContext outflow = (OutflowContext)ctx;
const NewtonianIdealGasContext gas = &outflow->gas;
const bool is_implicit = gas->is_implicit;

// Quadrature Point Loop
CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
// Setup
// -- Interp in
const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
State s = StateFromQ(gas, qi, state_var);
s.Y.pressure = outflow->pressure;
Expand All @@ -342,9 +336,11 @@ CEED_QFUNCTION_HELPER int PressureOutflow(void *ctx, CeedInt Q, const CeedScalar
for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j];

// Save values for Jacobian
StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur);
StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur);
} // End Quadrature Point Loop
if (is_implicit) {
StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur);
StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur);
}
}
return 0;
}

Expand Down
18 changes: 9 additions & 9 deletions examples/fluids/qfunctions/bc_slip.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,14 +13,11 @@
#include "riemann_solver.h"

CEED_QFUNCTION_HELPER int Slip(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) {
const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
const CeedScalar(*q_data_sur) = in[2];

CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
CeedScalar(*jac_data_sur) = out[1];

const NewtonianIdealGasContext newt_ctx = (const NewtonianIdealGasContext)ctx;
const CeedScalar zeros[6] = {0.};
const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
const CeedScalar(*q_data_sur) = in[2];
CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
CeedScalar(*jac_data_sur) = newt_ctx->is_implicit ? out[1] : NULL;

CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
Expand All @@ -42,8 +39,11 @@ CEED_QFUNCTION_HELPER int Slip(void *ctx, CeedInt Q, const CeedScalar *const *in
UnpackState_U(flux, Flux);
for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j];

StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur);
StoredValuesPack(Q, i, 5, 6, zeros, jac_data_sur); // Every output value must be set
if (newt_ctx->is_implicit) {
CeedScalar zeros[6] = {0.};
StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur);
StoredValuesPack(Q, i, 5, 6, zeros, jac_data_sur); // Every output value must be set
}
}
return 0;
}
Expand Down
Loading

0 comments on commit 780fa44

Please sign in to comment.