EB-FlowOverCylinder Case

Exec/RegTests/EB-FlowPastCylinder/CMakeLists.txt

@@ -0,0 +1,7 @@
+set(PELE_PHYSICS_EOS_MODEL GammaLaw)
+set(PELE_PHYSICS_CHEMISTRY_MODEL Null)
+set(PELE_PHYSICS_TRANSPORT_MODEL Constant)
+set(PELE_PHYSICS_ENABLE_SOOT OFF)
+set(PELE_PHYSICS_ENABLE_SPRAY OFF)
+set(PELE_PHYSICS_SPRAY_FUEL_NUM 0)
+include(BuildExeAndLib)

Exec/RegTests/EB-FlowPastCylinder/EB-FlowPastCylinder.inp

@@ -0,0 +1,77 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+# max_step = 100
+stop_time = 0.06
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic =  0  0  0
+geometry.coord_sys   =  0       # 0 => cart
+geometry.prob_lo     =  -2.   -2.  -2.
+geometry.prob_hi     =  10.  2.   2.
+amr.n_cell           =  192 64 16
+
+# >>>>>>>>>>>>>  BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<<
+# Interior, UserBC, Symmetry, SlipWall, NoSlipWall
+# >>>>>>>>>>>>>  BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<<
+
+pelec.lo_bc       =  "Hard" "SlipWall" "FOExtrap"
+pelec.hi_bc       =  "Hard" "SlipWall" "FOExtrap"
+
+# Problem setup
+prob.p = 1013250.
+prob.T = 298.
+prob.Re = 500.
+prob.Ma = 0.2
+prob.Pr = 0.7
+pelec.eb_boundary_T = 298.
+pelec.eb_isothermal = 1
+
+# WHICH PHYSICS
+pelec.do_hydro = 1
+pelec.do_mol = 0
+pelec.do_react = 0
+pelec.allow_negative_energy = 0
+pelec.diffuse_temp = 1
+pelec.diffuse_vel  = 1
+pelec.diffuse_spec = 0
+pelec.diffuse_enth = 1
+
+# TIME STEP CONTROL
+pelec.dt_cutoff      = 5.e-20  # level 0 timestep below which we halt
+pelec.cfl            = 0.7     # cfl number for hyperbolic system
+pelec.init_shrink    = 0.1    # scale back initial timestep
+pelec.change_max     = 1.05     # maximum increase in dt over successive steps
+
+# DIAGNOSTICS & VERBOSITY
+pelec.sum_interval   = 1       # timesteps between computing mass
+pelec.v              = 1       # verbosity in PeleC cpp files
+amr.v                = 1       # verbosity in Amr.cpp
+#amr.grid_log         = grdlog  # name of grid logging file
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 0       # maximum level number allowed
+amr.ref_ratio       = 2 2 2 2 # refinement ratio
+amr.regrid_int      = 5       # how often to regrid
+amr.blocking_factor = 16       # block factor in grid generation
+amr.max_grid_size   = 64
+
+# CHECKPOINT FILES
+amr.checkpoint_files_output = 0
+amr.check_file      = chk      # root name of checkpoint file
+amr.check_int       = 10000       # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_files_output = 1
+amr.plot_file       = plt
+amr.plot_int        = 5000
+amr.derive_plot_vars=ALL
+
+# EB
+ebd.boundary_grad_stencil_type = 0
+
+# for 2D (need a sphere):
+eb2.geom_type = "sphere"
+eb2.sphere_radius = 0.5
+eb2.sphere_center = 0 0 0
+eb2.sphere_has_fluid_inside = 0
+eb2.small_volfrac = 1.0e-4
+

Exec/RegTests/EB-FlowPastCylinder/GNUmakefile

@@ -0,0 +1,38 @@
+# AMReX
+DIM = 2
+COMP = gnu
+PRECISION = DOUBLE
+
+# Profiling
+PROFILE = FALSE
+TINY_PROFILE = FALSE
+COMM_PROFILE = FALSE
+TRACE_PROFILE = FALSE
+MEM_PROFILE = FALSE
+USE_GPROF = FALSE
+
+# Performance
+USE_MPI = TRUE
+USE_OMP = FALSE
+USE_CUDA = FALSE # invidia GPU
+USE_HIP = FALSE
+USE_SYCL = FALSE
+
+# Debugging
+DEBUG = FALSE
+FSANITIZER = FALSE
+THREAD_SANITIZER = FALSE
+
+# PeleC
+PELE_CVODE_FORCE_YCORDER = FALSE
+PELE_USE_MAGMA = FALSE
+PELE_COMPILE_AJACOBIAN = FALSE
+Eos_Model := GammaLaw
+Chemistry_Model := Null
+Transport_Model := Constant
+
+# GNU Make
+Bpack := ./Make.package
+Blocs := .
+PELE_HOME := ../../..
+include $(PELE_HOME)/Exec/Make.PeleC

Exec/RegTests/EB-FlowPastCylinder/Make.package

@@ -0,0 +1,3 @@
+CEXE_headers += prob.H
+CEXE_headers += prob_parm.H
+CEXE_sources += prob.cpp

Exec/RegTests/EB-FlowPastCylinder/README.md

@@ -0,0 +1,14 @@
+# 2D Hagen–Poiseuille flow over a cylinder
+
+Building on the [EB-C10 case](https://amrex-combustion.github.io/PeleC/ebverification/C10/README.html), this case prescribes Hagen-Poiseuille flow over a cylinder.
+
+
+## Short case description
+
+|                    | description                                         |
+|:-------------------|:----------------------------------------------------|
+| Problem definition | Hagen–Poiseuille flow over a cylinder               |
+| EB geometry        | embedded cylinder                                   |
+| EOS                | GammaLaw                                            |
+| Multi-level        | no                                                  |
+| Metric             | vortex shedding occurs                              |

Exec/RegTests/EB-FlowPastCylinder/prob.H

@@ -0,0 +1,144 @@
+#ifndef PROB_H
+#define PROB_H
+
+#include <AMReX_Geometry.H>
+#include <AMReX_FArrayBox.H>
+#include <AMReX_REAL.H>
+#include <AMReX_GpuMemory.H>
+
+#include "mechanism.H"
+
+#include "PeleC.H"
+#include "IndexDefines.H"
+#include "Constants.H"
+#include "PelePhysics.H"
+#include "Tagging.H"
+#include "ProblemSpecificFunctions.H"
+#include "prob_parm.H"
+
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+void
+pc_initdata(
+  int i,
+  int j,
+  int k,
+  amrex::Array4<amrex::Real> const& state,
+  amrex::GeometryData const& geomdata,
+  ProbParmDevice const& prob_parm)
+{
+  // Geometry
+  const amrex::Real* prob_lo = geomdata.ProbLo();
+  const amrex::Real* dx = geomdata.CellSize();
+
+  const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0];
+
+  const amrex::Real p = prob_parm.dpdx * x + prob_parm.p;
+  amrex::Real rho = 0.0, eint = 0.0;
+  amrex::Real massfrac[NUM_SPECIES] = {0.0};
+  for (int n = 0; n < NUM_SPECIES; n++) {
+    massfrac[n] = prob_parm.massfrac[n];
+  }
+  auto eos = pele::physics::PhysicsType::eos();
+  eos.PYT2RE(p, massfrac, prob_parm.T, rho, eint);
+
+  state(i, j, k, URHO) = rho;
+  state(i, j, k, UMX) = rho * prob_parm.uavg;
+  state(i, j, k, UMY) = 0.0;
+  state(i, j, k, UMZ) = 0.0;
+  state(i, j, k, UEINT) = rho * eint;
+  state(i, j, k, UEDEN) =
+    rho * (eint + 0.5 * (prob_parm.uavg * prob_parm.uavg));
+  state(i, j, k, UTEMP) = prob_parm.T;
+  for (int n = 0; n < NUM_SPECIES; n++) {
+    state(i, j, k, UFS + n) = rho * prob_parm.massfrac[n];
+  }
+}
+
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+void
+bcnormal(
+  const amrex::Real x[AMREX_SPACEDIM],
+  const amrex::Real s_inter[NVAR],
+  amrex::Real s_ext[NVAR],
+  const int idir,
+  const int sgn,
+  const amrex::Real /*time*/,
+  amrex::GeometryData const& geomdata,
+  ProbParmDevice const& prob_parm,
+  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& /*turb_fluc*/)
+{
+  if (idir == 0) {
+    amrex::Real rho = 0.0, u = 0.0, v = 0.0, w = 0.0, eint = 0.0, T = 0.0;
+    amrex::Real massfrac[NUM_SPECIES] = {0.0};
+    for (int n = 0; n < NUM_SPECIES; n++) {
+      massfrac[n] = prob_parm.massfrac[n];
+    }
+    auto eos = pele::physics::PhysicsType::eos();
+
+    if (sgn == 1) {
+      const amrex::Real p = prob_parm.dpdx * x[0] + prob_parm.p;
+      T = prob_parm.T;
+      eos.PYT2RE(p, massfrac, T, rho, eint);
+
+      u = s_inter[UMX] / s_inter[URHO];
+      v = s_inter[UMY] / s_inter[URHO];
+      w = s_inter[UMZ] / s_inter[URHO];
+
+    } else if (sgn == -1) {
+
+      // Following Blazek p 279, eq. 8.23
+
+      // Interior state (point d)
+      const amrex::Real* prob_hi = geomdata.ProbHi();
+      const amrex::Real* dx = geomdata.CellSize();
+      const amrex::Real xd = prob_hi[0] - 0.5 * dx[0];
+      const amrex::Real rho_inter = s_inter[URHO];
+      const amrex::Real u_inter = s_inter[UMX] / rho_inter;
+      const amrex::Real v_inter = s_inter[UMY] / rho_inter;
+      const amrex::Real w_inter = s_inter[UMZ] / rho_inter;
+      const amrex::Real T_inter = s_inter[UTEMP];
+      amrex::Real p_inter = 0.0, cs_inter = 0.0;
+      eos.RTY2P(rho_inter, T_inter, massfrac, p_inter);
+      eos.RTY2Cs(rho_inter, T_inter, massfrac, cs_inter);
+
+      // Boundary state (point b)
+      const amrex::Real xb = prob_hi[0];
+      const amrex::Real pb = prob_parm.dpdx * xb + prob_parm.p;
+      const amrex::Real rhob =
+        s_inter[URHO] + (pb - p_inter) / (cs_inter * cs_inter);
+      const amrex::Real ub = u_inter + (p_inter - pb) / (rho_inter * cs_inter);
+      const amrex::Real vb = v_inter;
+      const amrex::Real wb = w_inter;
+
+      // Ghost state (point a). Linear extrapolation from d and b
+      rho = (rhob - rho_inter) / (xb - xd) * (x[0] - xd) + rho_inter;
+      const amrex::Real p = (pb - p_inter) / (xb - xd) * (x[0] - xd) + p_inter;
+
+      eos.RYP2E(rho, massfrac, p, eint);
+      eos.EY2T(eint, massfrac, T);
+
+      u = (ub - u_inter) / (xb - xd) * (x[0] - xd) + u_inter;
+      v = (vb - v_inter) / (xb - xd) * (x[0] - xd) + v_inter;
+      w = (wb - w_inter) / (xb - xd) * (x[0] - xd) + w_inter;
+    }
+
+    s_ext[URHO] = rho;
+    s_ext[UMX] = rho * u;
+    s_ext[UMY] = rho * v;
+    s_ext[UMZ] = rho * w;
+    s_ext[UEINT] = rho * eint;
+    s_ext[UEDEN] = rho * (eint + 0.5 * (u * u + v * v + w * w));
+    s_ext[UTEMP] = T;
+    for (int n = 0; n < NUM_SPECIES; n++) {
+      s_ext[UFS + n] = rho * prob_parm.massfrac[n];
+    }
+  }
+}
+
+void pc_prob_close();
+
+using ProblemSpecificFunctions = DefaultProblemSpecificFunctions;
+
+#endif

Exec/RegTests/EB-FlowPastCylinder/prob.cpp

@@ -0,0 +1,100 @@
+#include "prob.H"
+
+void
+pc_prob_close()
+{
+}
+
+extern "C" {
+void
+amrex_probinit(
+  const int* /*init*/,
+  const int* /*name*/,
+  const int* /*namelen*/,
+  const amrex::Real* problo,
+  const amrex::Real* probhi)
+{
+  // Parse params
+  {
+    amrex::ParmParse pp("prob");
+    pp.query("p", PeleC::h_prob_parm_device->p);
+    pp.query("T", PeleC::h_prob_parm_device->T);
+    pp.query("Re", PeleC::h_prob_parm_device->Re);
+    pp.query("Ma", PeleC::h_prob_parm_device->Ma);
+    pp.query("Pr", PeleC::h_prob_parm_device->Pr);
+  }
+
+  amrex::Real L = (probhi[0] - problo[0]);
+  amrex::Real R = 0.5*(probhi[1] - problo[1]);
+  amrex::Real D = 2.0*R;
+
+  amrex::Real cp = 0.0;
+  amrex::Real cs = 0.0;
+  PeleC::h_prob_parm_device->massfrac[0] = 1.0;
+
+  auto eos = pele::physics::PhysicsType::eos();
+  eos.PYT2RE(
+    PeleC::h_prob_parm_device->p, PeleC::h_prob_parm_device->massfrac.begin(),
+    PeleC::h_prob_parm_device->T, PeleC::h_prob_parm_device->rho,
+    PeleC::h_prob_parm_device->eint);
+  eos.RTY2Cs(
+    PeleC::h_prob_parm_device->rho, PeleC::h_prob_parm_device->T,
+    PeleC::h_prob_parm_device->massfrac.begin(), cs);
+  eos.TY2Cp(
+    PeleC::h_prob_parm_device->T, PeleC::h_prob_parm_device->massfrac.begin(),
+    cp);
+
+  PeleC::h_prob_parm_device->umax = PeleC::h_prob_parm_device->Ma * cs;
+  PeleC::h_prob_parm_device->uavg = 0.5 * PeleC::h_prob_parm_device->umax;
+
+  auto& trans_parm = PeleC::trans_parms.host_parm();
+  trans_parm.const_bulk_viscosity = 0.0;
+  trans_parm.const_diffusivity = 0.0;
+  trans_parm.const_viscosity = PeleC::h_prob_parm_device->rho *
+                               PeleC::h_prob_parm_device->umax * D /
+                               PeleC::h_prob_parm_device->Re;
+  trans_parm.const_conductivity =
+    trans_parm.const_viscosity * cp / PeleC::h_prob_parm_device->Pr;
+  PeleC::trans_parms.sync_to_device();
+
+  PeleC::h_prob_parm_device->G =
+    PeleC::h_prob_parm_device->umax * 4 * trans_parm.const_viscosity /
+    (R * R);
+  PeleC::h_prob_parm_device->dpdx = -PeleC::h_prob_parm_device->G;
+
+  // Output IC
+  std::ofstream ofs("ic.txt", std::ofstream::out);
+  amrex::Print(ofs)
+    << "L, D, R, cs, rho, umax, p, T, gamma, mu, k, Re, Ma, Pr, dpdx, G"
+    << std::endl;
+  amrex::Print(ofs).SetPrecision(17)
+    << L << "," 
+    << D << ","
+    << R << ","
+    << cs << ","
+    << PeleC::h_prob_parm_device->rho << ","
+    << PeleC::h_prob_parm_device->umax << "," << PeleC::h_prob_parm_device->p
+    << "," << PeleC::h_prob_parm_device->T << "," << eos.gamma << ","
+    << trans_parm.const_viscosity << "," << trans_parm.const_conductivity << ","
+    << PeleC::h_prob_parm_device->Re << "," << PeleC::h_prob_parm_device->Ma
+    << "," << PeleC::h_prob_parm_device->Pr << ","
+    << PeleC::h_prob_parm_device->dpdx << "," << PeleC::h_prob_parm_device->G
+    << std::endl;
+  ofs.close();
+}
+}
+
+void
+PeleC::problem_post_timestep()
+{
+}
+
+void
+PeleC::problem_post_init()
+{
+}
+
+void
+PeleC::problem_post_restart()
+{
+}