Fix adCurvatures method

idaholab · May 24, 2023 · 6787636 · 6787636
1 parent 0c7e40b
commit 6787636
Show file tree

Hide file tree

Showing 4 changed files with 20 additions and 5 deletions.
diff --git a/framework/include/base/Assembly.h b/framework/include/base/Assembly.h
@@ -243,11 +243,7 @@ class Assembly
 
   const MooseArray<ADReal> & adJxWFace() const { return _ad_JxW_face; }
 
-  const MooseArray<ADReal> & adCurvatures() const
-  {
-    _calculate_curvatures = true;
-    return _ad_curvatures;
-  }
+  const MooseArray<ADReal> & adCurvatures() const;
 
   /**
    * Returns the reference to the coordinate transformation coefficients

diff --git a/framework/src/base/Assembly.C b/framework/src/base/Assembly.C
@@ -4736,6 +4736,19 @@ Assembly::processResidualsAndJacobian(const std::vector<ADReal> & residuals,
       cacheJacobian(row_indices[i], column_indices[j], element_matrix(i, j), matrix_tags);
 }
 
+const MooseArray<ADReal> &
+Assembly::adCurvatures() const
+{
+  _calculate_curvatures = true;
+  const Order helper_order = _mesh.hasSecondOrderElements() ? SECOND : FIRST;
+  const FEType helper_type(helper_order, LAGRANGE);
+  // Must prerequest the second derivatives. Sadly because there is only one _need_second_derivative
+  // map for both volumetric and face FE objects we must request both here
+  feSecondPhi<Real>(helper_type);
+  feSecondPhiFace<Real>(helper_type);
+  return _ad_curvatures;
+}
+
 template void coordTransformFactor<Point, Real>(const SubProblem & s,
                                                 SubdomainID sub_id,
                                                 const Point & point,

diff --git a/modules/navier_stokes/test/include/bcs/CrazyKCPlantFitsBoundary.h b/modules/navier_stokes/test/include/bcs/CrazyKCPlantFitsBoundary.h
@@ -47,6 +47,8 @@ class CrazyKCPlantFitsBoundary : public ADMaterial
   ADMaterialProperty<Real> & _surface_tension;
   ADMaterialProperty<RealVectorValue> & _grad_surface_tension;
   const MooseArray<ADPoint> & _ad_normals;
+  const MooseArray<ADReal> & _ad_curvatures;
+  ADMaterialProperty<RealVectorValue> & _surface_term_curvature;
   ADMaterialProperty<RealVectorValue> & _surface_term_gradient1;
   ADMaterialProperty<RealVectorValue> & _surface_term_gradient2;
 

diff --git a/modules/navier_stokes/test/src/bcs/CrazyKCPlantFitsBoundary.C b/modules/navier_stokes/test/src/bcs/CrazyKCPlantFitsBoundary.C
@@ -82,6 +82,8 @@ CrazyKCPlantFitsBoundary::CrazyKCPlantFitsBoundary(const InputParameters & param
     _grad_surface_tension(declareADProperty<RealVectorValue>(
         getParam<MaterialPropertyName>("grad_surface_tension_name"))),
     _ad_normals(_assembly.adNormals()),
+    _ad_curvatures(_assembly.adCurvatures()),
+    _surface_term_curvature(declareADProperty<RealVectorValue>("surface_term_curvature")),
     _surface_term_gradient1(declareADProperty<RealVectorValue>("surface_term_gradient1")),
     _surface_term_gradient2(declareADProperty<RealVectorValue>("surface_term_gradient2")),
     _length_units_per_meter(1. / std::pow(10, getParam<int>("length_unit_exponent"))),
@@ -115,6 +117,8 @@ CrazyKCPlantFitsBoundary::computeQpProperties()
   _grad_surface_tension[_qp] = _alpha * _mass_units_per_kilogram /
                                (_time_units_per_second * _time_units_per_second) /
                                _temperature_units_per_kelvin * _grad_temperature[_qp];
+  _surface_term_curvature[_qp] =
+      -2. * _ad_curvatures[_qp] * _surface_tension[_qp] * _ad_normals[_qp];
   _surface_term_gradient1[_qp] = -_grad_surface_tension[_qp];
   _surface_term_gradient2[_qp] = _ad_normals[_qp] * (_ad_normals[_qp] * _grad_surface_tension[_qp]);
 }