diff --git a/examples/fluids/README.md b/examples/fluids/README.md index 4a730a5d41..42a538adc1 100644 --- a/examples/fluids/README.md +++ b/examples/fluids/README.md @@ -441,8 +441,8 @@ For the 3D advection problem, the following additional command-line options are - `1E6` - `J` -* - `-bubble_type` - - `sphere` (3D) or `cylinder` (2D) +* - `-advection_ic_type` + - Initial condition type, from `sphere`, `cylinder`, `cosine_hill`, and `skew` - `sphere` - diff --git a/examples/fluids/advection.yaml b/examples/fluids/advection.yaml new file mode 100644 index 0000000000..42b210598d --- /dev/null +++ b/examples/fluids/advection.yaml @@ -0,0 +1,35 @@ +ceed: /cpu/self/memcheck/blocked +problem: advection +CtauS: .3 +stab: su +degree: 2 +dm_plex: + dim: 3 + box_faces: 5,5,1 + box_lower: 0,0,0 + box_upper: 125,125,250 + +bc_inflow: 1,2,3,4,5,6 + +units_kilogram: 1e-9 +rc: 100. +ksp: + atol: 1e-4 + rtol: 1e-3 + type: bcgs + +snes: + atol: 1e-3 + lag_jacobian: 100 + lag_jacobian_persists: + mf_operator: + +implicit: +ts: + dt: 1e-3 + type: alpha + max_steps: 10 + +dm_mat_preallocate_skip: 0 + +wind_type: rotation diff --git a/examples/fluids/navierstokes.c b/examples/fluids/navierstokes.c index 64875ca72e..8789196ef7 100644 --- a/examples/fluids/navierstokes.c +++ b/examples/fluids/navierstokes.c @@ -18,6 +18,8 @@ // ./navierstokes -ceed /cpu/self -problem density_current -degree 1 // ./navierstokes -ceed /gpu/cuda -problem advection -degree 1 // +//TESTARGS(name="Advection, skew") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/advection.yaml -ts_max_steps 0 -wind_type translation -wind_translation -0.5547002,0.83205029,0 -advection_ic_type skew -dm_plex_box_faces 2,1,1 -compare_final_state_atol 1e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-adv-skew.bin +//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 -pmat_pbdiagonal -pc_type vpbjacobi //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") -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 diff --git a/examples/fluids/navierstokes.h b/examples/fluids/navierstokes.h index 077ce564cc..39cd35f96a 100644 --- a/examples/fluids/navierstokes.h +++ b/examples/fluids/navierstokes.h @@ -55,8 +55,8 @@ static const char *const EulerTestTypes[] = {"isentropic_vortex", "test_1", // Advection - Wind types static const char *const WindTypes[] = {"rotation", "translation", "WindType", "WIND_", NULL}; -// Advection - Bubble Types -static const char *const BubbleTypes[] = {"sphere", "cylinder", "BubbleType", "BUBBLE_", NULL}; +// Advection - Initial Condition Types +static const char *const AdvectionICTypes[] = {"sphere", "cylinder", "cosine_hill", "skew", "AdvectionICType", "ADVECTIONIC_", NULL}; // Advection - Bubble Continuity Types static const char *const BubbleContinuityTypes[] = {"smooth", "back_sharp", "thick", "BubbleContinuityType", "BUBBLE_CONTINUITY_", NULL}; diff --git a/examples/fluids/problems/advection.c b/examples/fluids/problems/advection.c index 283bd1a3b8..f318216a5e 100644 --- a/examples/fluids/problems/advection.c +++ b/examples/fluids/problems/advection.c @@ -18,7 +18,7 @@ PetscErrorCode NS_ADVECTION(ProblemData *problem, DM dm, void *ctx, SimpleBC bc) { WindType wind_type; - BubbleType bubble_type; + AdvectionICType advectionic_type; BubbleContinuityType bubble_continuity_type; StabilizationType stab; SetupContextAdv setup_context; @@ -90,8 +90,8 @@ PetscErrorCode NS_ADVECTION(ProblemData *problem, DM dm, void *ctx, SimpleBC bc) PetscCall( PetscOptionsScalar("-strong_form", "Strong (1) or weak/integrated by parts (0) advection residual", NULL, strong_form, &strong_form, NULL)); PetscCall(PetscOptionsScalar("-E_wind", "Total energy of inflow wind", NULL, E_wind, &E_wind, NULL)); - PetscCall(PetscOptionsEnum("-bubble_type", "Sphere (3D) or cylinder (2D)", NULL, BubbleTypes, (PetscEnum)(bubble_type = BUBBLE_SPHERE), - (PetscEnum *)&bubble_type, NULL)); + PetscCall(PetscOptionsEnum("-advection_ic_type", "Initial condition for Advection problem", NULL, AdvectionICTypes, + (PetscEnum)(advectionic_type = ADVECTIONIC_BUBBLE_SPHERE), (PetscEnum *)&advectionic_type, NULL)); PetscCall(PetscOptionsEnum("-bubble_continuity", "Smooth, back_sharp, or thick", NULL, BubbleContinuityTypes, (PetscEnum)(bubble_continuity_type = BUBBLE_CONTINUITY_SMOOTH), (PetscEnum *)&bubble_continuity_type, NULL)); PetscCall(PetscOptionsEnum("-stab", "Stabilization method", NULL, StabilizationTypes, (PetscEnum)(stab = STAB_NONE), (PetscEnum *)&stab, NULL)); @@ -109,9 +109,10 @@ PetscErrorCode NS_ADVECTION(ProblemData *problem, DM dm, void *ctx, SimpleBC bc) if (wind_type == WIND_ROTATION && user_wind) { PetscCall(PetscPrintf(comm, "Warning! Use -wind_translation only with -wind_type translation\n")); } - if (wind_type == WIND_TRANSLATION && bubble_type == BUBBLE_CYLINDER && wind[2] != 0.) { + if (wind_type == WIND_TRANSLATION && advectionic_type == ADVECTIONIC_BUBBLE_CYLINDER && wind[2] != 0.) { wind[2] = 0; - PetscCall(PetscPrintf(comm, "Warning! Background wind in the z direction should be zero (-wind_translation x,x,0) with -bubble_type cylinder\n")); + PetscCall( + PetscPrintf(comm, "Warning! Background wind in the z direction should be zero (-wind_translation x,x,0) with -advection_ic_type cylinder\n")); } if (stab == STAB_NONE && CtauS != 0) { PetscCall(PetscPrintf(comm, "Warning! Use -CtauS only with -stab su or -stab supg\n")); @@ -154,7 +155,7 @@ PetscErrorCode NS_ADVECTION(ProblemData *problem, DM dm, void *ctx, SimpleBC bc) setup_context->wind[1] = wind[1]; setup_context->wind[2] = wind[2]; setup_context->wind_type = wind_type; - setup_context->bubble_type = bubble_type; + setup_context->initial_condition_type = advectionic_type; setup_context->bubble_continuity_type = bubble_continuity_type; setup_context->time = 0; @@ -194,12 +195,11 @@ PetscErrorCode PRINT_ADVECTION(User user, ProblemData *problem, AppCtx app_ctx) " Problem:\n" " Problem Name : %s\n" " Stabilization : %s\n" - " Bubble Type : %s (%" CeedInt_FMT "D)\n" + " Initial Condition Type : %s\n" " Bubble Continuity : %s\n" " Wind Type : %s\n", - app_ctx->problem_name, StabilizationTypes[advection_ctx->stabilization], BubbleTypes[setup_ctx->bubble_type], - setup_ctx->bubble_type == BUBBLE_SPHERE ? 3 : 2, BubbleContinuityTypes[setup_ctx->bubble_continuity_type], - WindTypes[setup_ctx->wind_type])); + app_ctx->problem_name, StabilizationTypes[advection_ctx->stabilization], AdvectionICTypes[setup_ctx->initial_condition_type], + BubbleContinuityTypes[setup_ctx->bubble_continuity_type], WindTypes[setup_ctx->wind_type])); if (setup_ctx->wind_type == WIND_TRANSLATION) { PetscCall(PetscPrintf(comm, " Background Wind : %f,%f,%f\n", setup_ctx->wind[0], setup_ctx->wind[1], setup_ctx->wind[2])); diff --git a/examples/fluids/qfunctions/advection.h b/examples/fluids/qfunctions/advection.h index 19ff253ae8..7979a0700a 100644 --- a/examples/fluids/qfunctions/advection.h +++ b/examples/fluids/qfunctions/advection.h @@ -14,8 +14,8 @@ #include #include -#include "../qfunctions/advection_types.h" -#include "../qfunctions/stabilization_types.h" +#include "advection_types.h" +#include "stabilization_types.h" #include "utils.h" typedef struct SetupContextAdv_ *SetupContextAdv; @@ -27,7 +27,7 @@ struct SetupContextAdv_ { CeedScalar wind[3]; CeedScalar time; WindType wind_type; - BubbleType bubble_type; + AdvectionICType initial_condition_type; BubbleContinuityType bubble_continuity_type; }; @@ -100,21 +100,27 @@ CEED_QFUNCTION_HELPER CeedInt Exact_Advection(CeedInt dim, CeedScalar time, cons // -- Energy CeedScalar r = 0.; - switch (context->bubble_type) { - case BUBBLE_SPHERE: { // (dim=3) + switch (context->initial_condition_type) { + case ADVECTIONIC_BUBBLE_SPHERE: // (dim=3) r = sqrt(Square(x - x0[0]) + Square(y - x0[1]) + Square(z - x0[2])); - } break; - case BUBBLE_CYLINDER: { // (dim=2) + break; + case ADVECTIONIC_BUBBLE_CYLINDER: // (dim=2) r = sqrt(Square(x - x0[0]) + Square(y - x0[1])); - } break; + break; + case ADVECTIONIC_COSINE_HILL: + r = sqrt(Square(x - center[0]) + Square(y - center[1])); + break; + case ADVECTIONIC_SKEW: + break; } // Initial Conditions + CeedScalar wind_scaling = 1.; switch (context->wind_type) { case WIND_ROTATION: q[0] = 1.; - q[1] = -(y - center[1]); - q[2] = (x - center[0]); + q[1] = -wind_scaling * (y - center[1]); + q[2] = wind_scaling * (x - center[0]); q[3] = 0; break; case WIND_TRANSLATION: @@ -125,20 +131,45 @@ CEED_QFUNCTION_HELPER CeedInt Exact_Advection(CeedInt dim, CeedScalar time, cons break; } - switch (context->bubble_continuity_type) { - // original continuous, smooth shape - case BUBBLE_CONTINUITY_SMOOTH: { - q[4] = r <= rc ? (1. - r / rc) : 0.; - } break; - // discontinuous, sharp back half shape - case BUBBLE_CONTINUITY_BACK_SHARP: { - q[4] = ((r <= rc) && (y < center[1])) ? (1. - r / rc) : 0.; + switch (context->initial_condition_type) { + case ADVECTIONIC_BUBBLE_SPHERE: + case ADVECTIONIC_BUBBLE_CYLINDER: + switch (context->bubble_continuity_type) { + // original continuous, smooth shape + case BUBBLE_CONTINUITY_SMOOTH: + q[4] = r <= rc ? (1. - r / rc) : 0.; + break; + // discontinuous, sharp back half shape + case BUBBLE_CONTINUITY_BACK_SHARP: + q[4] = ((r <= rc) && (y < center[1])) ? (1. - r / rc) : 0.; + break; + // attempt to define a finite thickness that will get resolved under grid refinement + case BUBBLE_CONTINUITY_THICK: + q[4] = ((r <= rc) && (y < center[1])) ? (1. - r / rc) * fmin(1.0, (center[1] - y) / 1.25) : 0.; + break; + } + break; + case ADVECTIONIC_COSINE_HILL: { + CeedScalar half_width = context->lx / 2; + q[4] = r > half_width ? 0. : cos(2 * M_PI * r / half_width + M_PI) + 1.; } break; - // attempt to define a finite thickness that will get resolved under grid refinement - case BUBBLE_CONTINUITY_THICK: { - q[4] = ((r <= rc) && (y < center[1])) ? (1. - r / rc) * fmin(1.0, (center[1] - y) / 1.25) : 0.; + case ADVECTIONIC_SKEW: { + CeedScalar skewed_barrier[3] = {wind[0], wind[1], 0}; + CeedScalar inflow_to_point[3] = {x - context->lx / 2, y, 0}; + CeedScalar cross_product[3] = {0}; + Cross3(skewed_barrier, inflow_to_point, cross_product); + + q[4] = cross_product[2] > 0 ? 0 : 1; + if ((x < 5 * CEED_EPSILON && wind[0] < 5 * CEED_EPSILON) || // outflow at -x boundary + (y < 5 * CEED_EPSILON && wind[1] < 5 * CEED_EPSILON) || // outflow at -y boundary + (x > context->lx - 5 * CEED_EPSILON && wind[0] > 5 * CEED_EPSILON) || // outflow at +x boundary + (y > context->ly - 5 * CEED_EPSILON && wind[1] > 5 * CEED_EPSILON) // outflow at +y boundary + ) { + q[4] = 0; + } } break; } + return 0; } diff --git a/examples/fluids/qfunctions/advection_types.h b/examples/fluids/qfunctions/advection_types.h index 4d62bea4fa..2d434cc2ff 100644 --- a/examples/fluids/qfunctions/advection_types.h +++ b/examples/fluids/qfunctions/advection_types.h @@ -16,11 +16,13 @@ typedef enum { WIND_TRANSLATION = 1, } WindType; -// Advection - Bubble Types +// Advection - Initial Condition Types typedef enum { - BUBBLE_SPHERE = 0, // dim=3 - BUBBLE_CYLINDER = 1, // dim=2 -} BubbleType; + ADVECTIONIC_BUBBLE_SPHERE = 0, // dim=3 + ADVECTIONIC_BUBBLE_CYLINDER = 1, // dim=2 + ADVECTIONIC_COSINE_HILL = 2, // dim=2 + ADVECTIONIC_SKEW = 3, +} AdvectionICType; // Advection - Bubble Continuity Types typedef enum { diff --git a/examples/fluids/qfunctions/shocktube.h b/examples/fluids/qfunctions/shocktube.h index 4064462227..8c10c43ff5 100644 --- a/examples/fluids/qfunctions/shocktube.h +++ b/examples/fluids/qfunctions/shocktube.h @@ -33,9 +33,6 @@ struct SetupContextShock_ { CeedScalar rho_high; CeedScalar P_low; CeedScalar rho_low; - int wind_type; // See WindType: 0=ROTATION, 1=TRANSLATION - int bubble_type; // See BubbleType: 0=SPHERE, 1=CYLINDER - int bubble_continuity_type; // See BubbleContinuityType: 0=SMOOTH, 1=BACK_SHARP 2=THICK }; typedef struct ShockTubeContext_ *ShockTubeContext; diff --git a/examples/fluids/tests-output/fluids-navierstokes-adv-rotation-cosine.bin b/examples/fluids/tests-output/fluids-navierstokes-adv-rotation-cosine.bin new file mode 100644 index 0000000000..295205a0b4 Binary files /dev/null and b/examples/fluids/tests-output/fluids-navierstokes-adv-rotation-cosine.bin differ diff --git a/examples/fluids/tests-output/fluids-navierstokes-adv-skew.bin b/examples/fluids/tests-output/fluids-navierstokes-adv-skew.bin new file mode 100644 index 0000000000..7a608cf4dc Binary files /dev/null and b/examples/fluids/tests-output/fluids-navierstokes-adv-skew.bin differ