Merge pull request CEED#1445 from CEED/jrwrigh/cosine_hill

fluids: Add cosine hill and skewed advection initial conditions

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`
   -
 

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

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

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};

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]));

examples/fluids/qfunctions/advection.h

@@ -14,8 +14,8 @@
 #include <ceed.h>
 #include <math.h>
 
-#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;
 }
 

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 {

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;