Merge pull request CEED#1557 from CEED/jrwrigh/diffusion_advection

fluids: Add diffusion to advection problem

examples/fluids/README.md

@@ -376,6 +376,11 @@ The following additional command-line options are available:
   - `1,0,0`
   -
 
+* - `-diffusion_coeff`
+  - Diffusion coefficient
+  - `0`
+  -
+
 * - `-E_wind`
   - Total energy of inflow wind when `-wind_type translation`
   - `1E6`

examples/fluids/index.md

@@ -505,15 +505,16 @@ The database may also be located on the same node as a MPI rank (collocated) or
 It's necessary to know how many ranks are associated with each collocated database, which is set by `-smartsim_collocated_database_num_ranks`.
 
 (problem-advection)=
-## Advection
+## Advection-Diffusion
 
 A simplified version of system {eq}`eq-ns`, only accounting for the transport of total energy, is given by
 
 $$
-\frac{\partial E}{\partial t} + \nabla \cdot (\bm{u} E ) = 0 \, ,
+\frac{\partial E}{\partial t} + \nabla \cdot (\bm{u} E ) - \kappa \nabla E = 0 \, ,
 $$ (eq-advection)
 
-with $\bm{u}$ the vector velocity field. In this particular test case, a blob of total energy (defined by a characteristic radius $r_c$) is transported by two different wind types.
+with $\bm{u}$ the vector velocity field and $\kappa$ the diffusion coefficient.
+In this particular test case, a blob of total energy (defined by a characteristic radius $r_c$) is transported by two different wind types.
 
 - **Rotation**
 

examples/fluids/navierstokes.c

@@ -20,7 +20,7 @@
 //
 //TESTARGS(name="Gaussian Wave, explicit, supg") -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 supg -state_var conservative -mass_ksp_type gmres -mass_pc_jacobi_type diagonal
 //TESTARGS(name="Advection 2D, rotation, explicit, supg, 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 supg -Ctaus 0.5 -mass_ksp_type gmres -mass_pc_type vpbjacobi -compare_final_state_atol 1e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-adv2d-rotation-explicit-stab-supg-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="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 -diffusion_coeff 5e-4 -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

examples/fluids/problems/advection.c

@@ -147,13 +147,14 @@ PetscErrorCode NS_ADVECTION(ProblemData *problem, DM dm, void *ctx, SimpleBC bc)
   // ------------------------------------------------------
   //             Create the libCEED context
   // ------------------------------------------------------
-  CeedScalar rc          = 1000.;  // m (Radius of bubble)
-  CeedScalar CtauS       = 0.;     // dimensionless
-  PetscBool  strong_form = PETSC_FALSE;
-  CeedScalar E_wind      = 1.e6;  // J
-  CeedScalar Ctau_a      = PetscPowScalarInt(user->app_ctx->degree, 2);
-  CeedScalar Ctau_t      = 0.;
-  PetscReal  wind[3]     = {1., 0, 0};  // m/s
+  CeedScalar rc              = 1000.;  // m (Radius of bubble)
+  CeedScalar CtauS           = 0.;     // dimensionless
+  PetscBool  strong_form     = PETSC_FALSE;
+  CeedScalar E_wind          = 1.e6;  // J
+  CeedScalar Ctau_a          = PetscPowScalarInt(user->app_ctx->degree, 2);
+  CeedScalar Ctau_t          = 0.;
+  PetscReal  wind[3]         = {1., 0, 0};  // m/s
+  CeedScalar diffusion_coeff = 0.;
   PetscReal  domain_min[3], domain_max[3], domain_size[3];
   PetscCall(DMGetBoundingBox(dm, domain_min, domain_max));
   for (PetscInt i = 0; i < problem->dim; i++) domain_size[i] = domain_max[i] - domain_min[i];
@@ -178,6 +179,7 @@ PetscErrorCode NS_ADVECTION(ProblemData *problem, DM dm, void *ctx, SimpleBC bc)
   PetscInt  n = problem->dim;
   PetscBool user_wind;
   PetscCall(PetscOptionsRealArray("-wind_translation", "Constant wind vector", NULL, wind, &n, &user_wind));
+  PetscCall(PetscOptionsScalar("-diffusion_coeff", "Diffusion coefficient", NULL, diffusion_coeff, &diffusion_coeff, NULL));
   PetscCall(PetscOptionsScalar("-CtauS", "Scale coefficient for tau (nondimensional)", NULL, CtauS, &CtauS, NULL));
   PetscCall(PetscOptionsBool("-strong_form", "Strong (true) or weak/integrated by parts (false) advection residual", NULL, strong_form, &strong_form,
                              NULL));
@@ -264,6 +266,7 @@ PetscErrorCode NS_ADVECTION(ProblemData *problem, DM dm, void *ctx, SimpleBC bc)
   advection_ctx->stabilization_tau = stab_tau;
   advection_ctx->Ctau_a            = Ctau_a;
   advection_ctx->Ctau_t            = Ctau_t;
+  advection_ctx->diffusion_coeff   = diffusion_coeff;
 
   PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &problem->ics.qfunction_context));
   PetscCallCeed(ceed,

examples/fluids/qfunctions/advection.h

@@ -326,6 +326,14 @@ CEED_QFUNCTION_HELPER void IFunction_AdvectionGeneric(void *ctx, CeedInt Q, cons
       for (CeedInt j = 0; j < dim; j++) grad_v[j][4][i] = -wdetJ * s.U.E_total * uX[j];
     }
 
+    {  // Diffusion
+      CeedScalar Fe[3], Fe_dXdx[3] = {0.};
+
+      for (CeedInt i = 0; i < dim; i++) Fe[i] = -context->diffusion_coeff * grad_s[i].U.E_total;
+      MatVecNM(dXdx, Fe, dim, dim, CEED_NOTRANSPOSE, Fe_dXdx);
+      for (CeedInt k = 0; k < dim; k++) grad_v[k][4][i] -= wdetJ * Fe_dXdx[k];
+    }
+
     const CeedScalar TauS = Tau(context, s, dXdx, dim);
     for (CeedInt j = 0; j < dim; j++) switch (context->stabilization) {
         case STAB_NONE:
@@ -453,6 +461,14 @@ CEED_QFUNCTION_HELPER void RHSFunction_AdvectionGeneric(void *ctx, CeedInt Q, co
       v[4][i] = 0.;
     }
 
+    {  // Diffusion
+      CeedScalar Fe[3], Fe_dXdx[3] = {0.};
+
+      for (CeedInt i = 0; i < dim; i++) Fe[i] = -context->diffusion_coeff * grad_s[i].U.E_total;
+      MatVecNM(dXdx, Fe, dim, dim, CEED_NOTRANSPOSE, Fe_dXdx);
+      for (CeedInt k = 0; k < dim; k++) grad_v[k][4][i] += wdetJ * Fe_dXdx[k];
+    }
+
     const CeedScalar TauS = Tau(context, s, dXdx, dim);
     for (CeedInt j = 0; j < dim; j++) switch (context->stabilization) {
         case STAB_NONE:

examples/fluids/qfunctions/advection_types.h

@@ -41,6 +41,7 @@ struct AdvectionContext_ {
   CeedScalar           Ctau_a;
   CeedScalar           Ctau_t;
   CeedScalar           dt;
+  CeedScalar           diffusion_coeff;
 };
 
 typedef struct SetupContextAdv_ *SetupContextAdv;