Merge pull request #1080 from CEED/ade/riemann-bcs

HLL Flux for Riemann BCs

Author: jedbrown

doc/sphinx/source/references.bib

@@ -202,3 +202,14 @@ @article{shurSTG2014
   doi = {10.1007/s10494-014-9534-8},
   langid = {english},
 }
+
+@inproceedings{mengaldoCompressibleBC2014,
+  title = {A Guide to the Implementation of Boundary Conditions in Compact High-Order Methods for Compressible Aerodynamics},
+  booktitle = {{{AIAA Aviation}} 2014},
+  author = {Mengaldo, Gianmarco and De Grazia, Daniele and Peiro, Joaquim and Farrington, Antony and Witherden, Freddie and Vincent, Peter and Sherwin, Spencer},
+  year = {2014},
+  month = jun,
+  publisher = {{AIAA}},
+  address = {{Atlanta}},
+  doi = {10.2514/6.2014-2923},
+}

examples/fluids/README.md

@@ -108,6 +108,10 @@ The following options are common among all problem types:
   - Use outflow boundary conditions on this list of faces
   -
 
+* - `-bc_freestream`
+  - Use freestream boundary conditions on this list of faces
+  -
+
 * - `-snes_view`
   - View PETSc `SNES` nonlinear solver configuration
   -
@@ -121,7 +125,7 @@ The following options are common among all problem types:
   -
 :::
 
-For the case of a square/cubic mesh, the list of face indices to be used with `-bc_wall`, `bc_inflow`, `bc_outflow` and/or `-bc_slip_x`, `-bc_slip_y`, and `-bc_slip_z` are:
+For the case of a square/cubic mesh, the list of face indices to be used with `-bc_wall`, `bc_inflow`, `bc_outflow`, `bc_freestream`  and/or `-bc_slip_x`, `-bc_slip_y`, and `-bc_slip_z` are:
 
 :::{list-table} 2D Face ID Labels
 :header-rows: 1
@@ -541,10 +545,63 @@ For the Density Current, Channel, and Blasius problems, the following common com
   - string
 :::
 
+#### Newtonian Wave
+
+The newtonian wave problem has the following command-line options in addition to the Newtonian Ideal Gas options:
+
+:::{list-table} Newtonian Wave Runtime Options
+:header-rows: 1
+
+* - Option
+  - Description
+  - Default value
+  - Unit
+
+* - `-velocity_freestream`
+  - Freestream velocity vector
+  - `0,0,0`
+  - `m/s`
+
+* - `-temperature_freestream`
+  - Freestream temperature
+  - `288`
+  - `K`
+
+* - `-pressure_freestream`
+  - Freestream pressure
+  - `1.01e5`
+  - `Pa`
+
+* - `-epicenter`
+  - Coordinates of center of perturbation
+  - `0,0,0`
+  - `m`
+
+* - `-amplitude`
+  - Amplitude of the perturbation
+  - `0.1`
+  -
+
+* - `-width`
+  - Width parameter of the perturbation
+  - `0.002`
+  - `m`
+
+:::
+
+This problem can be run with the `newtonianwave.yaml` file via:
+
+```
+./navierstokes -options_file newtonianwave.yaml
+```
+
+```{literalinclude} ../../../../../examples/fluids/newtonianwave.yaml
+:language: yaml
+```
 
 #### Density current
 
-The Density Current problem the following command-line options are available in addition to the Newtonian Ideal Gas options:
+The Density Current problem has the following command-line options in addition to the Newtonian Ideal Gas options:
 
 :::{list-table} Density Current Runtime Options
 :header-rows: 1
@@ -598,7 +655,7 @@ This problem can be run with:
 
 #### Channel flow
 
-The Channel problem the following command-line options are available in addition to the Newtonian Ideal Gas options:
+The Channel problem has the following command-line options in addition to the Newtonian Ideal Gas options:
 
 :::{list-table} Channel Runtime Options
 :header-rows: 1
@@ -640,7 +697,7 @@ This problem can be run with the `channel.yaml` file via:
 
 #### Blasius boundary layer
 
-The Blasius problem the following command-line options are available in addition to the Newtonian Ideal Gas options:
+The Blasius problem has the following command-line options in addition to the Newtonian Ideal Gas options:
 
 :::{list-table} Blasius Runtime Options
 :header-rows: 1

examples/fluids/index.md

@@ -385,7 +385,7 @@ $$
 \int_{\Omega} \nu_{SHOCK} \nabla \bm v \!:\! \nabla \bm q dV
 $$
 
-The shock capturing viscosity is implemented following the first formulation described in {cite} `tezduyar2007yzb`. The characteristic velocity $u_{cha}$ is taken to be the acoustic speed while the reference density $\rho_{ref}$ is just the local density. Shock capturing viscosity is defined by the following
+The shock capturing viscosity is implemented following the first formulation described in {cite}`tezduyar2007yzb`. The characteristic velocity $u_{cha}$ is taken to be the acoustic speed while the reference density $\rho_{ref}$ is just the local density. Shock capturing viscosity is defined by the following
 
 $$
 \nu_{SHOCK} = \tau_{SHOCK} u_{cha}^2
@@ -413,6 +413,23 @@ The constant $C_{YZB}$ is set to 0.1 for piecewise linear elements in the curren
 
 (problem-density-current)=
 
+## Newtonian Wave
+This test case is taken/inspired by that presented in {cite}`mengaldoCompressibleBC2014`. It is intended to test non-reflecting/Riemann boundary conditions. It's primarily intended for Euler equations, but has been implemented for the Navier-Stokes equations here for flexibility. 
+
+The problem has a perturbed initial condition and lets it evolve in time. The initial condition contains a Gaussian perturbation in the pressure field:
+
+$$
+\begin{aligned}
+\rho &= \rho_\infty\left(1+A\exp\left(\frac{-(\bar{x}^2 + \bar{y}^2)}{2\sigma^2}\right)\right) \\
+\bm{U} &= \bm U_\infty \\
+E &= \frac{p_\infty}{\gamma -1}\left(1+A\exp\left(\frac{-(\bar{x}^2 + \bar{y}^2)}{2\sigma^2}\right)\right) + \frac{\bm U_\infty \cdot \bm U_\infty}{2\rho_\infty},
+\end{aligned}
+$$
+
+where $A$ and $\sigma$ are the amplitude and width of the perturbation, respectively, and $(\bar{x}, \bar{y}) = (x-x_e, y-y_e)$ is the distance to the epicenter of the perturbation, $(x_e, y_e)$.
+
+Currently, a HLL (Harten, Lax, van Leer) Riemann solver is implemented under the freestream boundary condition. The formulation for the HLL solver is given in {cite}`toro2009`.
+
 ## Density Current
 
 For this test problem (from {cite}`straka1993numerical`), we solve the full Navier-Stokes equations {eq}`eq-ns`, for which a cold air bubble (of radius $r_c$) drops by convection in a neutrally stratified atmosphere.

examples/fluids/navierstokes.h

@@ -134,7 +134,8 @@ struct CeedData_private {
   CeedQFunction        qf_setup_vol, qf_ics, qf_rhs_vol, qf_ifunction_vol,
                        qf_setup_sur,
                        qf_apply_inflow, qf_apply_inflow_jacobian,
-                       qf_apply_outflow, qf_apply_outflow_jacobian;
+                       qf_apply_outflow, qf_apply_outflow_jacobian,
+                       qf_apply_freestream, qf_apply_freestream_jacobian;
   CeedBasis            basis_x, basis_xc, basis_q, basis_x_sur, basis_q_sur,
                        basis_xc_sur;
   CeedElemRestriction  elem_restr_x, elem_restr_q, elem_restr_qd_i;
@@ -182,8 +183,9 @@ struct SimpleBC_private {
             num_comps,
             num_slip[3], // Number of faces with slip BCs
             num_inflow,
-            num_outflow;
-  PetscInt  walls[16], slips[3][16], inflows[16], outflows[16];
+            num_outflow,
+            num_freestream;
+  PetscInt  walls[16], slips[3][16], inflows[16], outflows[16], freestreams[16];
   PetscBool user_bc;
 };
 
@@ -218,8 +220,8 @@ struct ProblemData_private {
   CeedInt           dim, q_data_size_vol, q_data_size_sur, jac_data_size_sur;
   CeedScalar        dm_scale;
   ProblemQFunctionSpec setup_vol, setup_sur, ics, apply_vol_rhs, apply_vol_ifunction,
-    apply_vol_ijacobian, apply_inflow, apply_outflow,
-    apply_inflow_jacobian, apply_outflow_jacobian;
+    apply_vol_ijacobian, apply_inflow, apply_outflow, apply_freestream,
+    apply_inflow_jacobian, apply_outflow_jacobian, apply_freestream_jacobian;
   bool              non_zero_time;
   PetscErrorCode    (*bc)(PetscInt, PetscReal, const PetscReal[], PetscInt,
                           PetscScalar[], void *);
@@ -235,6 +237,8 @@ extern int FreeContextPetsc(void *);
 // Set up problems
 // -----------------------------------------------------------------------------
 // Set up function for each problem
+extern PetscErrorCode NS_NEWTONIAN_WAVE(ProblemData *problem, DM dm,
+                                        void *ctx);
 extern PetscErrorCode NS_CHANNEL(ProblemData *problem, DM dm,
                                  void *ctx);
 extern PetscErrorCode NS_BLASIUS(ProblemData *problem, DM dm,
@@ -380,4 +384,5 @@ PetscErrorCode SetupStrongBC_Ceed(Ceed ceed, CeedData ceed_data, DM dm,
                                   User user, AppCtx app_ctx, ProblemData *problem,
                                   SimpleBC bc, CeedInt Q_sur, CeedInt q_data_size_sur);
 
+PetscErrorCode FreestreamBCSetup(ProblemData *problem, DM dm, void *ctx);
 #endif // libceed_fluids_examples_navier_stokes_h

examples/fluids/newtonianwave.yaml

@@ -0,0 +1,33 @@
+problem: 'newtonian_wave'
+mu: 0 # Effectively solving Euler momentum equations
+
+dm_plex_box_faces: 40,40,1
+dm_plex_box_upper: 1,1,0.1
+dm_plex_box_lower: 0,0,0
+dm_plex_dim: 3
+bc_freestream: 4,6,3,5
+bc_slip_z: 1,2
+
+velocity_freestream: 1,1,0
+temperature_freestream: 0.25
+pressure_freestream: 71.75
+
+epicenter: 0.33,0.75,0
+amplitude: 2
+width: 0.05
+
+ts:
+  adapt_type: 'none'
+  max_steps: 400
+  dt: 3.6e-4
+  type: alpha
+  alpha_radius: 0.5
+
+implicit: true
+stab: supg
+state_var: primitive
+
+snes_rtol: 1e-4
+ksp_rtol: 1e-2
+snes_lag_jacobian: 20
+snes_lag_jacobian_persists:

examples/fluids/problems/freestream_bc.c

@@ -0,0 +1,95 @@
+// Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors.
+// All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
+//
+// SPDX-License-Identifier: BSD-2-Clause
+//
+// This file is part of CEED:  http://github.com/ceed
+
+/// @file
+/// Utility functions for setting up Freestream boundary condition
+
+#include "../navierstokes.h"
+#include "../qfunctions/freestream_bc.h"
+
+PetscErrorCode FreestreamBCSetup(ProblemData *problem, DM dm, void *ctx) {
+  User      user    = *(User *)ctx;
+  MPI_Comm  comm    = PETSC_COMM_WORLD;
+  FreestreamContext freestream_ctx;
+  NewtonianIdealGasContext newtonian_ig_ctx;
+  CeedQFunctionContext freestream_context;
+  PetscFunctionBeginUser;
+
+  // *INDENT-OFF*
+  switch (user->phys->state_var) {
+  case STATEVAR_CONSERVATIVE:
+    problem->apply_freestream.qfunction              = Freestream_Conserv;
+    problem->apply_freestream.qfunction_loc          = Freestream_Conserv_loc;
+    problem->apply_freestream_jacobian.qfunction     = Freestream_Jacobian_Conserv;
+    problem->apply_freestream_jacobian.qfunction_loc = Freestream_Jacobian_Conserv_loc;
+  case STATEVAR_PRIMITIVE:
+    problem->apply_freestream.qfunction              = Freestream_Prim;
+    problem->apply_freestream.qfunction_loc          = Freestream_Prim_loc;
+    problem->apply_freestream_jacobian.qfunction     = Freestream_Jacobian_Prim;
+    problem->apply_freestream_jacobian.qfunction_loc = Freestream_Jacobian_Prim_loc;
+  }
+  // *INDENT-ON*
+
+  // -- Option Defaults
+  CeedScalar U_inf[3] = {0.};   // m/s
+  CeedScalar T_inf    = 288.;   // K
+  CeedScalar P_inf    = 1.01e5; // Pa
+
+  PetscOptionsBegin(comm, NULL, "Options for Freestream boundary condition",
+                    NULL);
+  PetscInt narray=3;
+  PetscCall(PetscOptionsScalarArray("-velocity_freestream",
+                                    "Velocity at freestream condition",
+                                    NULL, U_inf, &narray, NULL));
+  PetscCheck(narray == 3, comm, PETSC_ERR_ARG_SIZ,
+             "-velocity_freestream should recieve array of size 3, instead recieved size %"
+             PetscInt_FMT".", narray);
+
+  PetscCall(PetscOptionsScalar("-temperature_freestream",
+                               "Temperature at freestream condition",
+                               NULL, T_inf, &T_inf, NULL));
+  PetscCall(PetscOptionsScalar("-pressure_freestream",
+                               "Pressure at freestream condition",
+                               NULL, P_inf, &P_inf, NULL));
+  PetscOptionsEnd();
+
+  PetscScalar meter  = user->units->meter;
+  PetscScalar second = user->units->second;
+  PetscScalar Kelvin = user->units->Kelvin;
+  PetscScalar Pascal = user->units->Pascal;
+
+  T_inf *= Kelvin;
+  P_inf *= Pascal;
+  for (int i=0; i<3; i++) U_inf[i] *= meter / second;
+
+  CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context,
+                              CEED_MEM_HOST, &newtonian_ig_ctx);
+  State S_infty;
+  {
+    CeedScalar Y[5] = {P_inf, U_inf[0], U_inf[1], U_inf[2], T_inf};
+    CeedScalar x[3] = {0.};
+    S_infty = StateFromY(newtonian_ig_ctx, Y, x);
+  }
+
+  // -- Set freestream_ctx struct values
+  PetscCall(PetscCalloc1(1, &freestream_ctx));
+  freestream_ctx->newtonian_ctx = *newtonian_ig_ctx;
+  freestream_ctx->S_infty       = S_infty;
+
+  CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context,
+                                  &newtonian_ig_ctx);
+
+  CeedQFunctionContextCreate(user->ceed, &freestream_context);
+  CeedQFunctionContextSetData(freestream_context, CEED_MEM_HOST, CEED_USE_POINTER,
+                              sizeof(*freestream_ctx), freestream_ctx);
+  CeedQFunctionContextSetDataDestroy(freestream_context, CEED_MEM_HOST,
+                                     FreeContextPetsc);
+  problem->apply_freestream.qfunction_context = freestream_context;
+  problem->apply_freestream_jacobian.qfunction_context = freestream_context;
+
+  PetscFunctionReturn(0);
+};