[WIP] Dome natural convection

.gitignore

@@ -83,6 +83,7 @@ temp_print_trace.*
 *.net
 .depend
 *.png
+*.avi
 *.svg
 *.gif
 *.tif

framework/include/linearfvbcs/LinearFVAdvectionDiffusionFunctorDirichletBC.h

@@ -43,4 +43,6 @@ class LinearFVAdvectionDiffusionFunctorDirichletBC : public LinearFVAdvectionDif
 protected:
   /// The functor for this BC (can be variable, function, etc)
   const Moose::Functor<Real> & _functor;
+
+  const enum NormalComponent : int { X = 0, Y, Z, NONE } _normal_component;
 };

framework/include/timeintegrators/BDF2.h

@@ -30,6 +30,10 @@ class BDF2 : public TimeIntegrator
                                 ADReal & ad_u_dotdot) const override;
   virtual void postResidual(NumericVector<Number> & residual) override;
   virtual bool overridesSolve() const override { return false; }
+  virtual Real timeDerivativeRHSContribution(const dof_id_type dof_id,
+                                             const std::vector<Real> & factors) const override;
+  virtual Real timeDerivativeMatrixContribution(const Real factor) const override;
+  virtual unsigned int numStatesRequired() const override { return 2; }
 
 protected:
   /**

framework/src/kokkos/systems/KokkosAuxiliarySystem.K

@@ -29,7 +29,7 @@ AuxiliarySystem::addKokkosKernel(const std::string & kernel_name,
 void
 AuxiliarySystem::kokkosCompute(ExecFlagType type)
 {
-  TIME_SECTION("computeKokkosAuxKernel", 1);
+  TIME_SECTION("computeKokkosAuxKernel", 3);
 
   if (!_kokkos_elemental_aux_storage[type].size() && !_kokkos_nodal_aux_storage[type].size())
     return;

framework/src/linearfvbcs/LinearFVAdvectionDiffusionFunctorDirichletBC.C

@@ -20,19 +20,29 @@ LinearFVAdvectionDiffusionFunctorDirichletBC::validParams()
       "finite volume system and whose face values are determined using a functor. This kernel is "
       "only designed to work with advection-diffusion problems.");
   params.addRequiredParam<MooseFunctorName>("functor", "The functor for this boundary condition.");
+  MooseEnum normal_component("x y z none", "none");
+  params.addParam<MooseEnum>("normal_component",
+                             normal_component,
+                             "What normal component to multiply the Dirichlet value by (no "
+                             "multiplication occurs if 'none' is selected)");
   return params;
 }
 
 LinearFVAdvectionDiffusionFunctorDirichletBC::LinearFVAdvectionDiffusionFunctorDirichletBC(
     const InputParameters & parameters)
-  : LinearFVAdvectionDiffusionBC(parameters), _functor(getFunctor<Real>("functor"))
+  : LinearFVAdvectionDiffusionBC(parameters),
+    _functor(getFunctor<Real>("functor")),
+    _normal_component((getParam<MooseEnum>("normal_component").getEnum<NormalComponent>()))
 {
 }
 
 Real
 LinearFVAdvectionDiffusionFunctorDirichletBC::computeBoundaryValue() const
 {
-  return _functor(singleSidedFaceArg(_current_face_info), determineState());
+  auto diri_value = _functor(singleSidedFaceArg(_current_face_info), determineState());
+  if (_normal_component != NONE)
+    diri_value *= _current_face_info->normal()(_normal_component);
+  return diri_value;
 }
 
 Real
@@ -42,9 +52,7 @@ LinearFVAdvectionDiffusionFunctorDirichletBC::computeBoundaryNormalGradient() co
                                         ? _current_face_info->elemPtr()
                                         : _current_face_info->neighborPtr());
   const Real distance = computeCellToFaceDistance();
-  return (_functor(singleSidedFaceArg(_current_face_info), determineState()) -
-          raw_value(_var(elem_arg, determineState()))) /
-         distance;
+  return (computeBoundaryValue() - raw_value(_var(elem_arg, determineState()))) / distance;
 }
 
 Real
@@ -59,7 +67,7 @@ Real
 LinearFVAdvectionDiffusionFunctorDirichletBC::computeBoundaryValueRHSContribution() const
 {
   // Fetch the boundary value from the provided functor.
-  return _functor(singleSidedFaceArg(_current_face_info), determineState());
+  return computeBoundaryValue();
 }
 
 Real
@@ -75,6 +83,5 @@ LinearFVAdvectionDiffusionFunctorDirichletBC::computeBoundaryGradientRHSContribu
 {
   // The boundary term from the central difference approximation of the
   // normal gradient.
-  return _functor(singleSidedFaceArg(_current_face_info), determineState()) /
-         computeCellToFaceDistance();
+  return computeBoundaryValue() / computeCellToFaceDistance();
 }

framework/src/systems/LinearSystem.C

@@ -367,7 +367,12 @@ LinearSystem::containsTimeKernel()
 }
 
 void
-LinearSystem::compute(ExecFlagType)
+LinearSystem::compute(const ExecFlagType type)
 {
-  // Linear systems have their own time derivative computation machinery
+  // - Linear system assembly is associated with EXEC_NONLINEAR
+  // - Avoid division by 0 dt
+  if (type == EXEC_NONLINEAR && _fe_problem.dt() > 0.)
+    for (auto & ti : _time_integrators)
+      // Do things like compute integration weights
+      ti->preStep();
 }

framework/src/systems/SolverSystem.C

@@ -150,10 +150,11 @@ SolverSystem::compute(const ExecFlagType type)
       compute_tds = true;
   }
 
+  // avoid division by dt which might be zero.
   if (compute_tds && _fe_problem.dt() > 0.)
     for (auto & ti : _time_integrators)
     {
-      // avoid division by dt which might be zero.
+      // Do things like compute integration weights
       ti->preStep();
       ti->computeTimeDerivatives();
     }

framework/src/timeintegrators/BDF2.C

@@ -85,3 +85,26 @@ BDF2::duDotDuCoeff() const
   else
     return _weight[0];
 }
+
+Real
+BDF2::timeDerivativeRHSContribution(dof_id_type dof_id, const std::vector<Real> & factors) const
+{
+  mooseAssert(factors.size() == numStatesRequired(),
+              "Either too many or too few states are given!");
+
+  if (_t_step == 1)
+    return factors[0] * _solution_old(dof_id) / _dt;
+  else
+    return -(_weight[1] * factors[0] * _solution_old(dof_id) +
+             _weight[2] * factors[1] * _solution_older(dof_id)) /
+           _dt;
+}
+
+Real
+BDF2::timeDerivativeMatrixContribution(const Real factor) const
+{
+  if (_t_step == 1)
+    return factor / _dt;
+  else
+    return factor * _weight[0] / _dt;
+}

modules/navier_stokes/include/auxkernels/LinearWCNSFVSmagorinskyTurbulentViscosityAux.h

@@ -0,0 +1,44 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "AuxKernel.h"
+
+class INSFVVelocityVariable;
+
+/*
+ * Computes the value of the eddy viscosity for the mixing length model.
+ */
+class LinearWCNSFVSmagorinskyTurbulentViscosityAux : public AuxKernel
+{
+public:
+  static InputParameters validParams();
+
+  LinearWCNSFVSmagorinskyTurbulentViscosityAux(const InputParameters & parameters);
+
+protected:
+  virtual Real computeValue() override;
+
+  /// the dimension of the simulation
+  const unsigned int _dim;
+
+  /// x-velocity
+  MooseLinearVariableFVReal * const _u_var;
+  /// y-velocity
+  MooseLinearVariableFVReal * const _v_var;
+  /// z-velocity
+  MooseLinearVariableFVReal * const _w_var;
+
+  /// Density
+  const Moose::Functor<Real> & _rho;
+
+  /// Smagorinsky constant
+  const Real _Cs;
+};

modules/navier_stokes/src/auxkernels/LinearWCNSFVSmagorinskyTurbulentViscosityAux.C

@@ -0,0 +1,90 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "LinearWCNSFVSmagorinskyTurbulentViscosityAux.h"
+#include "NS.h"
+
+registerMooseObject("NavierStokesApp", LinearWCNSFVSmagorinskyTurbulentViscosityAux);
+
+InputParameters
+LinearWCNSFVSmagorinskyTurbulentViscosityAux::validParams()
+{
+  InputParameters params = AuxKernel::validParams();
+  params.addClassDescription("Computes the turbulent viscosity for the mixing length model.");
+  params.addRequiredParam<SolverVariableName>("u", "The velocity in the x direction.");
+  params.addParam<SolverVariableName>("v", "The velocity in the y direction.");
+  params.addParam<SolverVariableName>("w", "The velocity in the z direction.");
+  params.addParam<Real>("smagorinsky_constant", 0.18, "Value of Smagorinsky's constant to use");
+  params.addRequiredParam<MooseFunctorName>(NS::density, "Density functor");
+  return params;
+}
+
+LinearWCNSFVSmagorinskyTurbulentViscosityAux::LinearWCNSFVSmagorinskyTurbulentViscosityAux(
+    const InputParameters & params)
+  : AuxKernel(params),
+    _dim(_subproblem.mesh().dimension()),
+    _u_var(dynamic_cast<MooseLinearVariableFVReal *>(
+        &_subproblem.getVariable(_tid, getParam<SolverVariableName>("u")))),
+    _v_var(params.isParamValid("v")
+               ? dynamic_cast<MooseLinearVariableFVReal *>(
+                     &_subproblem.getVariable(_tid, getParam<SolverVariableName>("v")))
+               : nullptr),
+    _w_var(params.isParamValid("w")
+               ? dynamic_cast<MooseLinearVariableFVReal *>(
+                     &_subproblem.getVariable(_tid, getParam<SolverVariableName>("w")))
+               : nullptr),
+    _rho(getFunctor<Real>(NS::density)),
+    _Cs(getParam<Real>("smagorinsky_constant"))
+{
+  if (!_u_var)
+    paramError("u", "the u velocity must be a MooseLinearVariableFVReal.");
+
+  if (_dim >= 2 && !_v_var)
+    paramError("v",
+               "In two or more dimensions, the v velocity must be supplied and it must be a "
+               "MooseLinearVariableFVReal.");
+
+  if (_dim >= 3 && !_w_var)
+    paramError("w",
+               "In three-dimensions, the w velocity must be supplied and it must be a "
+               "MooseLinearVariableFVReal.");
+
+  _u_var->computeCellGradients();
+  if (_v_var)
+    _v_var->computeCellGradients();
+  if (_w_var)
+    _w_var->computeCellGradients();
+}
+
+Real
+LinearWCNSFVSmagorinskyTurbulentViscosityAux::computeValue()
+{
+  const auto & elem_info = _mesh.elemInfo(_current_elem->id());
+
+  const auto grad_u = _u_var->gradSln(elem_info);
+  Real symmetric_strain_tensor_norm = 2.0 * Utility::pow<2>(grad_u(0));
+  if (_dim >= 2)
+  {
+    const auto grad_v = _v_var->gradSln(elem_info);
+    symmetric_strain_tensor_norm +=
+        2.0 * Utility::pow<2>(grad_v(1)) + Utility::pow<2>(grad_v(0) + grad_u(1));
+    if (_dim >= 3)
+    {
+      const auto grad_w = _w_var->gradSln(elem_info);
+      symmetric_strain_tensor_norm += 2.0 * Utility::pow<2>(grad_w(2)) +
+                                      Utility::pow<2>(grad_u(2) + grad_w(0)) +
+                                      Utility::pow<2>(grad_v(2) + grad_w(1));
+    }
+  }
+
+  symmetric_strain_tensor_norm = std::sqrt(symmetric_strain_tensor_norm);
+
+  return _rho(makeElemArg(_current_elem), determineState()) * symmetric_strain_tensor_norm *
+         libMesh::Utility::pow<2>(_Cs * _current_elem->hmax());
+}

modules/navier_stokes/src/userobjects/RhieChowMassFlux.C

@@ -511,8 +511,8 @@ RhieChowMassFlux::populateCouplingFunctors(
             {
               force_kernel->setCurrentElemInfo(&elem_info);
               face_hbya(dim_i) -= force_kernel->computeRightHandSideContribution() *
-                                  ainv_reader[dim_i](elem_dof) /
-                                  elem_info.volume(); // zero-term expansion
+                                  ainv_reader[dim_i](elem_dof) / elem_info.volume() /
+                                  elem_info.coordFactor(); // zero-term expansion
             }
           face_hbya(dim_i) *= boundary_normal_multiplier;
         }