Diffuse interface support for more IC's and IC user objects idaholab#…

…13331

framework/include/ics/BoundingBoxIC.h

@@ -61,5 +61,7 @@ class BoundingBoxIC : public InitialCondition
 
   /// The Point object constructed from the x2, y2, z2 components for the bottom left BB corner
   const Point _top_right;
+
+  const Real _int_width;
 };
 

framework/src/ics/BoundingBoxIC.C

@@ -9,6 +9,7 @@
 
 #include "BoundingBoxIC.h"
 #include "libmesh/point.h"
+#include "libmesh/utility.h"
 
 registerMooseObject("MooseApp", BoundingBoxIC);
 
@@ -28,6 +29,9 @@ validParams<BoundingBoxIC>()
   params.addParam<Real>("inside", 0.0, "The value of the variable inside the box");
   params.addParam<Real>("outside", 0.0, "The value of the variable outside the box");
 
+  params.addParam<Real>(
+      "int_width", 0.0, "The width of the diffuse interface. Set to 0 for sharp interface.");
+
   params.addClassDescription("BoundingBoxIC allows setting the initial condition of a value inside "
                              "and outside of a specified box. The box is aligned with the x, y, z "
                              "axes");
@@ -46,16 +50,30 @@ BoundingBoxIC::BoundingBoxIC(const InputParameters & parameters)
     _inside(getParam<Real>("inside")),
     _outside(getParam<Real>("outside")),
     _bottom_left(_x1, _y1, _z1),
-    _top_right(_x2, _y2, _z2)
+    _top_right(_x2, _y2, _z2),
+    _int_width(getParam<Real>("int_width"))
 {
 }
 
 Real
 BoundingBoxIC::value(const Point & p)
 {
-  for (unsigned int i = 0; i < LIBMESH_DIM; i++)
-    if (p(i) < _bottom_left(i) || p(i) > _top_right(i))
-      return _outside;
+  if (_int_width <= 0.0)
+  {
+    for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
+      if (p(i) < _bottom_left(i) || p(i) > _top_right(i))
+        return _outside;
+
+    return _inside;
+  }
+  else
+  {
+    Real f_in = 1.0;
+    for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
+      if (_bottom_left(i) != _top_right(i))
+        f_in *= 0.5 * (std::tanh(2.0 * (p(i) - _bottom_left(i)) / _int_width) -
+                       std::tanh(libMesh::pi * (p(i) - _top_right(i)) / _int_width));
 
-  return _inside;
+    return _outside + (_inside - _outside) * f_in;
+  }
 }

modules/phase_field/include/userobjects/PolycrystalCircles.h

@@ -48,9 +48,13 @@ class PolycrystalCircles : public PolycrystalUserObjectBase
     R
   }; // Names of columns in text file.
   const bool _columnar_3D;
+  Real _int_width;         // Width of diffuse interface
   unsigned int _grain_num; // Number of crystal grains to create
 
   std::vector<Point> _centerpoints; // x,y,z coordinates of circle centers
   std::vector<Real> _radii;         // Radius for each circular grain created
+
+private:
+  Real computeDiffuseInterface(const Point & p, const unsigned int & i) const;
 };
 

modules/phase_field/include/userobjects/PolycrystalVoronoi.h

@@ -37,6 +37,7 @@ class PolycrystalVoronoi : public PolycrystalUserObjectBase
   const bool _columnar_3D;
 
   const unsigned int _rand_seed;
+  const Real _int_width;
 
   Point _bottom_left;
   Point _top_right;
@@ -45,5 +46,16 @@ class PolycrystalVoronoi : public PolycrystalUserObjectBase
   std::vector<Point> _centerpoints;
 
   const FileName _file_name;
+
+private:
+  Real computeDiffuseInterface(const Point & point,
+                               const unsigned int & gr_index,
+                               const std::vector<unsigned int> & grain_ids) const;
+  Point findNormalVector(const Point & point, const Point & p1, const Point & p2) const;
+  Point findCenterPoint(const Point & point, const Point & p1, const Point & p2) const;
+  Real findLinePoint(const Point & point,
+                     const Point & N,
+                     const Point & cntr,
+                     const unsigned int dim) const;
 };
 

modules/phase_field/src/ics/SmoothCircleBaseIC.C

@@ -137,8 +137,7 @@ SmoothCircleBaseIC::computeCircleValue(const Point & p, const Point & center, co
     }
 
     case ProfileType::TANH:
-      return (_invalue - _outvalue) * 0.5 *
-                 (std::tanh(libMesh::pi * (radius - dist) / _int_width) + 1.0) +
+      return (_invalue - _outvalue) * 0.5 * (std::tanh(2.0 * (radius - dist) / _int_width) + 1.0) +
              _outvalue;
 
     default:

modules/phase_field/src/userobjects/PolycrystalCircles.C

@@ -11,6 +11,7 @@
 #include "MooseMesh.h"
 #include "MooseVariable.h"
 
+#include "libmesh/utility.h"
 #include <fstream>
 
 registerMooseObject("PhaseFieldApp", PolycrystalCircles);
@@ -34,13 +35,15 @@ validParams<PolycrystalCircles>()
   params.addParam<std::vector<Real>>("z_positions", "z coordinate for each circle center");
   params.addParam<std::vector<Real>>("radii", "The radius for each circle");
   params.addParam<FileName>("file_name", "File containing circle centers and radii");
+  params.addParam<Real>("int_width", 0, "Width of diffuse interface");
 
   return params;
 }
 
 PolycrystalCircles::PolycrystalCircles(const InputParameters & parameters)
   : PolycrystalUserObjectBase(parameters),
     _columnar_3D(getParam<bool>("columnar_3D")),
+    _int_width(getParam<Real>("int_width")),
     _grain_num(0)
 {
 }
@@ -65,7 +68,7 @@ PolycrystalCircles::getGrainsBasedOnPoint(const Point & point,
     else
       distance = _mesh.minPeriodicDistance(_vars[0]->number(), _centerpoints[i], point);
 
-    if (distance < _radii[i])
+    if (distance < _radii[i] + _int_width)
       grains.push_back(i);
   }
 }
@@ -84,7 +87,7 @@ PolycrystalCircles::getVariableValue(unsigned int op_index, const Point & p) con
       break;
     }
 
-  return active_grain_on_op != invalid_id ? 1.0 : 0.0;
+  return active_grain_on_op != invalid_id ? computeDiffuseInterface(p, active_grain_on_op) : 0.0;
 }
 
 void
@@ -170,3 +173,23 @@ PolycrystalCircles::precomputeGrainStructure()
     }
   }
 }
+
+Real
+PolycrystalCircles::computeDiffuseInterface(const Point & p, const unsigned int & i) const
+{
+  if (_int_width == 0)
+    return 1.0;
+
+  Real d = 0;
+
+  if (_columnar_3D)
+  {
+    Real d_x = (p(0) - _centerpoints[i](0)) * (p(0) - _centerpoints[i](0));
+    Real d_y = (p(1) - _centerpoints[i](1)) * (p(1) - _centerpoints[i](1));
+    d = std::sqrt(d_x + d_y);
+  }
+  else
+    d = _mesh.minPeriodicDistance(_vars[0]->number(), _centerpoints[i], p);
+
+  return 0.5 * (1 - std::tanh(2.0 * (d - _radii[i]) / _int_width));
+}

modules/phase_field/src/userobjects/PolycrystalVoronoi.C

@@ -34,6 +34,7 @@ validParams<PolycrystalVoronoi>()
       "",
       "File containing grain centroids, if file_name is provided, the centroids "
       "from the file will be used.");
+  params.addParam<Real>("int_width", 0.0, "Width of diffuse interfaces");
   return params;
 }
 
@@ -42,6 +43,7 @@ PolycrystalVoronoi::PolycrystalVoronoi(const InputParameters & parameters)
     _grain_num(getParam<unsigned int>("grain_num")),
     _columnar_3D(getParam<bool>("columnar_3D")),
     _rand_seed(getParam<unsigned int>("rand_seed")),
+    _int_width(getParam<Real>("int_width")),
     _file_name(getParam<FileName>("file_name"))
 {
   if (_file_name == "" && _grain_num == 0)
@@ -54,7 +56,7 @@ PolycrystalVoronoi::PolycrystalVoronoi(const InputParameters & parameters)
 void
 PolycrystalVoronoi::getGrainsBasedOnPoint(const Point & point,
                                           std::vector<unsigned int> & grains) const
-{
+{ /*
   auto n_grains = _centerpoints.size();
   auto min_distance = _range.norm();
   auto min_index = n_grains;
@@ -74,7 +76,39 @@ PolycrystalVoronoi::getGrainsBasedOnPoint(const Point & point,
   mooseAssert(min_index < n_grains, "Couldn't find closest Voronoi cell");
 
   grains.resize(1);
-  grains[0] = min_index;
+  grains[0] = min_index; */
+
+  grains.resize(0);
+  Real d_min = _range.norm();
+  Real distance;
+  auto n_grains = _centerpoints.size();
+  auto min_index = n_grains;
+
+  for (MooseIndex(_centerpoints) grain = 0; grain < n_grains; ++grain)
+  {
+    distance = _mesh.minPeriodicDistance(_vars[0]->number(), _centerpoints[grain], point);
+    if (distance < d_min)
+    {
+      d_min = distance;
+      min_index = grain;
+    }
+  }
+  mooseAssert(min_index < n_grains, "Couldn't find closest Voronoi cell");
+
+  Point current_grain = _centerpoints[min_index];
+  grains.push_back(min_index);
+  for (MooseIndex(_centerpoints) grain = 0; grain < n_grains; ++grain)
+  {
+    if (grain != min_index)
+    {
+      Point next_grain = _centerpoints[grain];
+      Point N = findNormalVector(point, current_grain, next_grain);
+      Point cntr = findCenterPoint(point, current_grain, next_grain);
+      distance = N * (cntr - point);
+      if (distance < 0.6 * _int_width)
+        grains.push_back(grain);
+    }
+  }
 }
 
 Real
@@ -92,7 +126,11 @@ PolycrystalVoronoi::getVariableValue(unsigned int op_index, const Point & p) con
       break;
     }
 
-  return active_grain_on_op != invalid_id ? 1.0 : 0.0;
+  Real profile_val = 0.0;
+  if (active_grain_on_op != invalid_id)
+    profile_val = computeDiffuseInterface(p, active_grain_on_op, grain_ids);
+
+  return profile_val;
 }
 
 void
@@ -151,3 +189,60 @@ PolycrystalVoronoi::precomputeGrainStructure()
     }
   }
 }
+
+Real
+PolycrystalVoronoi::computeDiffuseInterface(const Point & point,
+                                            const unsigned int & gr_index,
+                                            const std::vector<unsigned int> & grain_ids) const
+{
+  Real val = 1.0;
+  Point current_grain = _centerpoints[gr_index];
+  for (auto i : grain_ids)
+    if (i != gr_index)
+    {
+      Point next_grain = _centerpoints[i];
+      Point N = findNormalVector(point, current_grain, next_grain);
+      Point cntr = findCenterPoint(point, current_grain, next_grain);
+      for (unsigned int dim = 0; dim < 3; ++dim)
+        if (N(dim) != 0.0)
+        {
+          Real L = findLinePoint(point, N, cntr, dim);
+          val *= 0.5 * (1.0 - std::tanh(2.0 * (point(dim) - L) * N(dim) / _int_width));
+          break;
+        }
+    }
+  return val;
+}
+
+Point
+PolycrystalVoronoi::findNormalVector(const Point & point, const Point & p1, const Point & p2) const
+{
+  Point pa = point + _mesh.minPeriodicVector(_vars[0]->number(), point, p1);
+  Point pb = point + _mesh.minPeriodicVector(_vars[0]->number(), point, p2);
+  Point N = pb - pa;
+  return N / N.norm();
+}
+
+Point
+PolycrystalVoronoi::findCenterPoint(const Point & point, const Point & p1, const Point & p2) const
+{
+  Point pa = point + _mesh.minPeriodicVector(_vars[0]->number(), point, p1);
+  Point pb = point + _mesh.minPeriodicVector(_vars[0]->number(), point, p2);
+  return 0.5 * (pa + pb);
+}
+
+Real
+PolycrystalVoronoi::findLinePoint(const Point & point,
+                                  const Point & N,
+                                  const Point & cntr,
+                                  const unsigned int dim) const
+{
+  Real Lsum;
+  if (dim == 0)
+    Lsum = N(1) * (point(1) - cntr(1)) + N(2) * (point(2) - cntr(2));
+  else if (dim == 1)
+    Lsum = N(0) * (point(0) - cntr(0)) + N(2) * (point(2) - cntr(2));
+  else
+    Lsum = N(0) * (point(0) - cntr(0)) + N(1) * (point(1) - cntr(1));
+  return cntr(dim) - Lsum / N(dim);
+}

modules/phase_field/test/tests/DeformedGrain/DeformedGrain.i

@@ -40,6 +40,7 @@
     type = PolycrystalVoronoi
     rand_seed = 81
     coloring_algorithm = bt
+    int_width = 0
   [../]
   [./grain_tracker]
     type = GrainTracker

modules/phase_field/test/tests/MultiPhase/switchingfunctionmultiphasematerial.i

@@ -10,7 +10,6 @@
 []
 
 [GlobalParams]
-  int_width = 3.0
   block = 0
 []
 
@@ -40,6 +39,7 @@
     outside = 0.0
     type = BoundingBoxIC
     variable = eta2
+    int_width = 0
   [../]
   [./IC_eta3]
     x1 = 15
@@ -50,6 +50,7 @@
     outside = 0.0
     type = BoundingBoxIC
     variable = eta3
+    int_width = 0
   [../]
   [./IC_eta4]
     x1 = 0
@@ -60,6 +61,7 @@
     outside = 0.0
     type = BoundingBoxIC
     variable = eta0
+    int_width = 0
   [../]
   [./IC_c]
     x1 = 15
@@ -69,6 +71,7 @@
     variable = c
     invalue = 1.0
     type = SmoothCircleIC
+    int_width = 3.0
   [../]
   [./IC_eta1]
     x1 = 15
@@ -78,6 +81,7 @@
     variable = eta1
     invalue = 1.0
     type = SmoothCircleIC
+    int_width = 3.0
   [../]
 []
 

modules/phase_field/test/tests/initial_conditions/PolycrystalVoronoiVoidIC_moregrains.i

@@ -46,6 +46,7 @@
 [UserObjects]
   [./voronoi]
     type = PolycrystalVoronoi
+    int_width = 0
   [../]
 []
 

modules/phase_field/test/tests/initial_conditions/PolycrystalVoronoiVoidIC_notperiodic.i

@@ -50,6 +50,7 @@
   [./voronoi]
     type = PolycrystalVoronoi
     rand_seed = 12444
+    int_width = 0
   [../]
 []
 

modules/phase_field/test/tests/initial_conditions/PolycrystalVoronoiVoidIC_periodic.i

@@ -50,6 +50,7 @@
   [./voronoi]
     type = PolycrystalVoronoi
     rand_seed = 10
+    int_width = 0
   [../]
 []