rotating reference frames capability

CMake/BuildPeleExe.cmake

@@ -81,6 +81,7 @@ function(build_pele_exe pele_exe_name pele_physics_lib_name)
        ${SRC_DIR}/Utilities.H
        ${SRC_DIR}/Utilities.cpp
        ${SRC_DIR}/WENO.H
+       ${SRC_DIR}/RotSource.cpp
   )
 
   if(PELE_PHYSICS_ENABLE_SPRAY)

Docs/sphinx/geometry/EB.rst

@@ -127,7 +127,7 @@ the end of each time step. The re-redistribution is performed every time the ref
 presented in `Pember et al. <https://www.sciencedirect.com/science/article/pii/S0021999185711655>`_. A forthcoming paper will describe the methodology for this procedure when using state redistribution.
 
 
-Date Structures and utility functions
+Data Structures and utility functions
 -------------------------------------
 
 Several structures exist to store geometry dependent information. These are populated on creation of a new AMRLevel and stored in the PeleC object so that they are available for computation. These facilitate accessing the EB data. The datatypes are:
@@ -205,3 +205,90 @@ Problem specific inflow conditions on an EB
 It is possible for the user to define problem specific conditions on an EB surface. This is done by defining an ``problem_eb_state`` function and then including `pelec.eb_problem_state = 1` in the input file. An example of this is found in the ``EB-InflowBC`` case.
 
 .. warning:: This is a beta feature. Currently this will only affect the calculation of the hydrodynamic fluxes so this works best for advection dominated EB conditions.
+
+Rotational frames
+------------------------------
+
+.. warning:: PeleC has limited support for simulations using rotating frames using the single-reference-frame (SRF) 
+   method only in 3D and needs the MOL scheme when using this feature.
+
+
+Single-reference-frame implementation
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+PeleC currently supports a non-inertial frame based implementation 
+for rotating bodies, which are relevant for turbomachinery applications. 
+The source terms for momentum and energy equations
+are applied in the non-inertial reference frame for cylindrically symmetric
+embedded boundaries. The total energy in the non-inertial frame is 
+modified with the inclusion of rotatational kinetic energy:
+
+.. math::
+   E = e + \frac{u^2+v^2+w^2}{2} - \frac{\omega^2 r^2}{2} 
+
+where E is the total specific energy, e is specific internal energy,
+u,v,w are the velocity components, :math:`\omega` is the angular velocity, and 
+r is the distance from the rotational axis.
+The momentum equation includes additional source terms for 
+Coriolis, centrifugal, and rotational acceleration:
+
+.. math::
+
+        \frac{\partial}{\partial t} (\rho \mathbf{v})+ \mathbf{\nabla} \cdot (\rho V V)=-\mathbf{\nabla}P + \mathbf{\nabla} \cdot \bar{\bar{\tau}}
+        -\rho (2 \mathbf{\omega} \times \mathbf{V} + \mathbf{\omega} \times (\mathbf{\omega} \times r) + \dot{\mathbf{\omega}} \times r)
+
+For further details, refer to Appendix A.3 
+*Navier-Stokes Equations in Rotating Frame of Reference* 
+in textbook *Computational fluid dynamics: principles and applications* by J. Blazek. 
+
+Single-reference-frame usage
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Most often in turbomachinery applications, there is rotor and 
+a stator, which will be represented in PeleC as embedded boundaries.
+The rotor will most often include blades and other components that are 
+not rotational symmetric like a cylindrical stator surface. This capability
+will solve the governing equations in the non inertial frame where the 
+rotor boundary would then be a no-slip wall while the stator, which is rotationally 
+symmetric, will be a slip wall. Therefore, the user in this case
+has to differentiate where the stator is using ``pelec.rf_rad`` parameter that 
+says anything beyond this radius will be a stator.
+The rotational frame capability can be turned on by using these inputs:
+
+::
+
+   pelec.do_rf=1  #turn on rotational frames
+   pelec.rf_omega=500.0 #rotational speed (rad/s)
+   pelec.rf_axis=2      #axis of rotation (0-x,1-y,2-z)
+   pelec.rf_axis_x=0.0  #x location of axis 
+   pelec.rf_axis_y=0.0  #y location of axis
+   pelec.rf_axis_z=0.0  #z location of axis
+   pelec.rf_rad=20.0    #non-zero EB velocity beyond this radius
+   pelec.do_mol = 1     #rotational frames only implemented with MOL
+
+The velocities in the plot files will be in the rotating frame.
+But, one can turn on ``amr.derive_plot_vars`` to add new variables 
+that transform velocity back to inertial frame. These variable names 
+include ``x_velocity_if``, ``y_velocity_if`` , ``z_velocity_if`` and 
+``magvel_if``.
+
+Some example test cases
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+We have tested this SRF implementation using a Taylor-Couette example (see ``Exec/EB_TaylorCouette``)
+for which we obtain solutions matching the analytic solution as shown in
+Figure :ref:`rotational-frame-verification <rotframe>` (a).
+
+We have also done a comparison with the incompressible SRF solver in OpenFOAM for a
+4 bladed rotor case (see ``Exec/EB_Rotor4Blade``). 
+The solutions agree reasonably well, as shown in Figure :ref:`rotational-frame-verification <rotframe>` (b) 
+with minor differences due to the compressible formulation in PeleC. 
+
+.. _rotframe:
+
+.. figure:: rotframetests.png
+   :alt: EB Cell
+   :width: 400
+
+   \(a) solution to Taylor-Couette flow compared to analytic solution and (b) 
+   comparison of velocity magntiude for a 4 blade rotor case with incompressible 
+   SRF solver from OpenFOAM. 

Exec/RegTests/EB_Rotor4Blade/GNUmakefile

@@ -0,0 +1,38 @@
+# AMReX
+DIM = 3
+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 = FALSE
+USE_OMP = FALSE
+USE_CUDA = FALSE
+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_Rotor4Blade/Make.package

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

Exec/RegTests/EB_Rotor4Blade/prob.H

@@ -0,0 +1,171 @@
+#ifndef PROB_H
+#define PROB_H
+
+#include <AMReX_Print.H>
+#include <AMReX_ParmParse.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 "EOS.H"
+
+#include "Tagging.H"
+#include "ProblemSpecificFunctions.H"
+#include "prob_parm.H"
+#include "EB.H"
+#include <AMReX_EB2_IF_Rotation.H>
+#include "Utilities.H"
+#include "Geometry.H"
+
+class EBConcCylinders : public pele::pelec::Geometry::Register<EBConcCylinders>
+{
+public:
+  static std::string identifier() { return "conc_cylinders"; }
+
+  void
+  build(const amrex::Geometry& geom, const int max_coarsening_level) override;
+};
+
+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();
+
+  amrex::Real u = 0.0;
+  amrex::Real v = 0.0;
+  amrex::Real w = 0.0;
+  const amrex::Real p = prob_parm.p0;
+  const amrex::Real T = prob_parm.T0;
+  amrex::Real massfrac[NUM_SPECIES] = {1.0};
+  amrex::Real rho, eint;
+  auto eos = pele::physics::PhysicsType::eos();
+  eos.PYT2RE(p, massfrac, T, rho, eint);
+
+  if (prob_parm.do_rf) {
+    amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> axis_loc = {
+      prob_parm.rf_axis_x, prob_parm.rf_axis_y, prob_parm.rf_axis_z};
+    amrex::RealVect r = get_rotaxis_vec(i, j, k, axis_loc, geomdata);
+    amrex::RealVect omgavec(0.0, 0.0, 0.0);
+    omgavec[prob_parm.rf_axis] = prob_parm.rf_omega;
+    amrex::RealVect w_cross_r = omgavec.crossProduct(r);
+    u = -w_cross_r[0];
+    v = -w_cross_r[1];
+    w = -w_cross_r[2];
+  }
+
+  state(i, j, k, URHO) = rho;
+  state(i, j, k, UMX) = rho * u;
+  state(i, j, k, UMY) = rho * v;
+  state(i, j, k, UMZ) = rho * w;
+  state(i, j, k, UEINT) = rho * eint;
+  state(i, j, k, UEDEN) = rho * (eint + 0.5 * (u * u + v * v + w * w));
+  state(i, j, k, UTEMP) = T;
+  for (int n = 0; n < NUM_SPECIES; n++) {
+    state(i, j, k, UFS + n) = rho * massfrac[n];
+  }
+
+  if (prob_parm.do_rf) {
+    amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> axis_loc = {
+      prob_parm.rf_axis_x, prob_parm.rf_axis_y, prob_parm.rf_axis_z};
+    amrex::Real rotenrg = get_rot_energy(
+      i, j, k, prob_parm.rf_omega, prob_parm.rf_axis, axis_loc, geomdata);
+    state(i, j, k, UEDEN) -= rho * rotenrg;
+  }
+}
+
+AMREX_GPU_DEVICE
+AMREX_FORCE_INLINE
+void
+bcnormal(
+  const amrex::Real* /*x[AMREX_SPACEDIM]*/,
+  const amrex::Real* /*s_int[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*/)
+{
+}
+
+void pc_prob_close();
+
+using ProblemSpecificFunctions = DefaultProblemSpecificFunctions;
+
+/*using ProblemTags = EmptyProbTagStruct;
+
+  struct MyProbDeriveStruct
+  {
+  static amrex::Box the_same_box(const amrex::Box& b) { return b; }
+
+  static void
+  add(amrex::DeriveList& derive_lst, const amrex::DescriptorList& desc_lst)
+  {
+// Add derives as follows and define the derive function below:
+derive_lst.add(
+"velx_if", amrex::IndexType::TheCellType(), 1, pc_velx_if,
+the_same_box);
+derive_lst.addComponent("velx_if", desc_lst, State_Type, 0, NVAR);
+
+derive_lst.add(
+"vely_if", amrex::IndexType::TheCellType(), 1, pc_vely_if,
+the_same_box);
+derive_lst.addComponent("vely_if", desc_lst, State_Type, 0, NVAR);
+
+derive_lst.add(
+"velz_if", amrex::IndexType::TheCellType(), 1, pc_velz_if,
+the_same_box);
+derive_lst.addComponent("velz_if", desc_lst, State_Type, 0, NVAR);
+}
+
+static void pc_rhoerror(
+const amrex::Box& bx,
+amrex::FArrayBox& derfab,
+int dcomp,
+int ncomp,
+const amrex::FArrayBox& datfab,
+const amrex::Geometry& geomdata,
+amrex::Real time,
+const int* bcrec,
+int level)
+{
+auto const dat = datfab.array();
+auto arr = derfab.array();
+const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> prob_lo =
+geomdata.ProbLoArray();
+const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx =
+geomdata.CellSizeArray();
+
+ProbParmDevice const* prob_parm = PeleC::d_prob_parm_device;
+
+amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0];
+const amrex::Real rho_exact =
+prob_parm->rho +
+prob_parm->alpha * std::exp(
+-((x - prob_parm->cs * time) / prob_parm->sigma) *
+((x - prob_parm->cs * time) / prob_parm->sigma));
+arr(i, j, k, 0) = dat(i, j, k, URHO) - rho_exact;
+});
+}
+}; */
+
+#endif