Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "SmoothCircleBaseIC.h"
      11             : #include "MooseMesh.h"
      12             : #include "MooseVariable.h"
      13             : #include "FEProblemBase.h"
      14             : 
      15             : #include "libmesh/utility.h"
      16             : 
      17             : InputParameters
      18        5468 : SmoothCircleBaseIC::validParams()
      19             : {
      20        5468 :   InputParameters params = InitialCondition::validParams();
      21       10936 :   params.addRequiredParam<Real>("invalue", "The variable value inside the circle");
      22       10936 :   params.addRequiredParam<Real>("outvalue", "The variable value outside the circle");
      23       10936 :   params.addParam<Real>(
      24       10936 :       "int_width", 0.0, "The interfacial width of the void surface.  Defaults to sharp interface");
      25       10936 :   params.addParam<bool>("3D_spheres", true, "in 3D, whether the objects are spheres or columns");
      26       10936 :   params.addParam<bool>("zero_gradient",
      27       10936 :                         false,
      28             :                         "Set the gradient DOFs to zero. This can avoid "
      29             :                         "numerical problems with higher order shape "
      30             :                         "functions and overlapping circles.");
      31       10936 :   params.addParam<unsigned int>("rand_seed", 12345, "Seed value for the random number generator");
      32       10936 :   MooseEnum profileType("COS TANH", "COS");
      33       10936 :   params.addParam<MooseEnum>(
      34             :       "profile", profileType, "Functional dependence for the interface profile");
      35        5468 :   return params;
      36        5468 : }
      37             : 
      38        2849 : SmoothCircleBaseIC::SmoothCircleBaseIC(const InputParameters & parameters)
      39             :   : InitialCondition(parameters),
      40        2849 :     _mesh(_fe_problem.mesh()),
      41        2849 :     _invalue(parameters.get<Real>("invalue")),
      42        2849 :     _outvalue(parameters.get<Real>("outvalue")),
      43        2849 :     _int_width(parameters.get<Real>("int_width")),
      44        2849 :     _3D_spheres(parameters.get<bool>("3D_spheres")),
      45        2849 :     _zero_gradient(parameters.get<bool>("zero_gradient")),
      46        2942 :     _num_dim(_3D_spheres ? 3 : 2),
      47        8547 :     _profile(getParam<MooseEnum>("profile").getEnum<ProfileType>())
      48             : {
      49        5698 :   _random.seed(_tid, getParam<unsigned int>("rand_seed"));
      50             : 
      51        2849 :   if (_int_width <= 0.0 && _profile == ProfileType::TANH)
      52           0 :     paramError("int_width",
      53             :                "Interface width has to be strictly positive for the hyperbolic tangent profile");
      54        2849 : }
      55             : 
      56             : void
      57        2158 : SmoothCircleBaseIC::initialSetup()
      58             : {
      59             :   // Compute radii and centers and initialize vector sizes
      60        2158 :   computeCircleRadii();
      61        2158 :   computeCircleCenters();
      62             : 
      63        2158 :   if (_centers.size() != _radii.size())
      64           0 :     mooseError("_center and _radii vectors are not the same size in the Circle IC");
      65             : 
      66        2158 :   if (_centers.size() < 1)
      67           0 :     mooseError("_center and _radii were not initialized in the Circle IC");
      68        2158 : }
      69             : 
      70             : Real
      71    13543496 : SmoothCircleBaseIC::value(const Point & p)
      72             : {
      73    13543496 :   Real value = _outvalue;
      74             :   Real val2 = 0.0;
      75             : 
      76   331826366 :   for (unsigned int circ = 0; circ < _centers.size() && value != _invalue; ++circ)
      77             :   {
      78   318282870 :     val2 = computeCircleValue(p, _centers[circ], _radii[circ]);
      79   318282870 :     if ((val2 > value && _invalue > _outvalue) || (val2 < value && _outvalue > _invalue))
      80             :       value = val2;
      81             :   }
      82             : 
      83    13543496 :   return value;
      84             : }
      85             : 
      86             : RealGradient
      87      228588 : SmoothCircleBaseIC::gradient(const Point & p)
      88             : {
      89      228588 :   if (_zero_gradient)
      90             :     return 0.0;
      91             : 
      92             :   RealGradient gradient = 0.0;
      93      228588 :   Real value = _outvalue;
      94             :   Real val2 = 0.0;
      95             : 
      96      502176 :   for (unsigned int circ = 0; circ < _centers.size(); ++circ)
      97             :   {
      98      273588 :     val2 = computeCircleValue(p, _centers[circ], _radii[circ]);
      99      273588 :     if ((val2 > value && _invalue > _outvalue) || (val2 < value && _outvalue > _invalue))
     100             :     {
     101             :       value = val2;
     102       77754 :       gradient = computeCircleGradient(p, _centers[circ], _radii[circ]);
     103             :     }
     104             :   }
     105             : 
     106      228588 :   return gradient;
     107             : }
     108             : 
     109             : Real
     110   318496458 : SmoothCircleBaseIC::computeCircleValue(const Point & p, const Point & center, const Real & radius)
     111             : {
     112   318496458 :   Point l_center = center;
     113   318496458 :   Point l_p = p;
     114   318496458 :   if (!_3D_spheres) // Create 3D cylinders instead of spheres
     115             :   {
     116      216288 :     l_p(2) = 0.0;
     117      216288 :     l_center(2) = 0.0;
     118             :   }
     119             :   // Compute the distance between the current point and the center
     120   318496458 :   Real dist = _mesh.minPeriodicDistance(_var.number(), l_p, l_center);
     121             : 
     122   318496458 :   switch (_profile)
     123             :   {
     124   318446858 :     case ProfileType::COS:
     125             :     {
     126             :       // Return value
     127   318446858 :       Real value = _outvalue; // Outside circle
     128             : 
     129   318446858 :       if (dist <= radius - _int_width / 2.0) // Inside circle
     130     4012601 :         value = _invalue;
     131   314434257 :       else if (dist < radius + _int_width / 2.0) // Smooth interface
     132             :       {
     133     1600904 :         Real int_pos = (dist - radius + _int_width / 2.0) / _int_width;
     134     1600904 :         value = _outvalue + (_invalue - _outvalue) * (1.0 + std::cos(int_pos * libMesh::pi)) / 2.0;
     135             :       }
     136             :       return value;
     137             :     }
     138             : 
     139       49600 :     case ProfileType::TANH:
     140       49600 :       return (_invalue - _outvalue) * 0.5 * (std::tanh(2.0 * (radius - dist) / _int_width) + 1.0) +
     141       49600 :              _outvalue;
     142             : 
     143           0 :     default:
     144           0 :       mooseError("Internal error.");
     145             :   }
     146             : }
     147             : 
     148             : RealGradient
     149       77754 : SmoothCircleBaseIC::computeCircleGradient(const Point & p,
     150             :                                           const Point & center,
     151             :                                           const Real & radius)
     152             : {
     153       77754 :   Point l_center = center;
     154       77754 :   Point l_p = p;
     155       77754 :   if (!_3D_spheres) // Create 3D cylinders instead of spheres
     156             :   {
     157           0 :     l_p(2) = 0.0;
     158           0 :     l_center(2) = 0.0;
     159             :   }
     160             :   // Compute the distance between the current point and the center
     161       77754 :   Real dist = _mesh.minPeriodicDistance(_var.number(), l_p, l_center);
     162             : 
     163             :   // early return if we are probing the center of the circle
     164       77754 :   if (dist == 0.0)
     165             :     return 0.0;
     166             : 
     167             :   Real DvalueDr = 0.0;
     168       77672 :   switch (_profile)
     169             :   {
     170       47672 :     case ProfileType::COS:
     171       47672 :       if (dist < radius + _int_width / 2.0 && dist > radius - _int_width / 2.0)
     172             :       {
     173       18276 :         const Real int_pos = (dist - radius + _int_width / 2.0) / _int_width;
     174       18276 :         const Real Dint_posDr = 1.0 / _int_width;
     175       18276 :         DvalueDr = Dint_posDr * (_invalue - _outvalue) *
     176       18276 :                    (-std::sin(int_pos * libMesh::pi) * libMesh::pi) / 2.0;
     177             :       }
     178             :       break;
     179             : 
     180       30000 :     case ProfileType::TANH:
     181       30000 :       DvalueDr = -(_invalue - _outvalue) * 0.5 / _int_width * libMesh::pi *
     182       30000 :                  (1.0 - Utility::pow<2>(std::tanh(4.0 * (radius - dist) / _int_width)));
     183       30000 :       break;
     184             : 
     185           0 :     default:
     186           0 :       mooseError("Internal error.");
     187             :   }
     188             : 
     189       77672 :   return _mesh.minPeriodicVector(_var.number(), center, p) * (DvalueDr / dist);
     190             : }