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 "FourierNoise.h"
      11             : #include "MooseRandom.h"
      12             : #include "FEProblemBase.h"
      13             : #include "libmesh/utility.h"
      14             : 
      15             : registerMooseObject("PhaseFieldApp", FourierNoise);
      16             : 
      17             : InputParameters
      18          23 : FourierNoise::validParams()
      19             : {
      20          23 :   InputParameters params = Function::validParams();
      21          23 :   params.addClassDescription("Generate noise from a fourier series");
      22          46 :   params.addRequiredParam<Real>("lambda",
      23             :                                 "Wavelength cut off (set to about twice the interfacial width)");
      24          46 :   params.addParam<unsigned int>(
      25             :       "num_terms",
      26             :       "Number of random fourier series terms (this will result in non-periodic noise). Omit this "
      27             :       "parameter to obtain a periodic noise distribution.");
      28          46 :   params.addParam<unsigned int>("seed", 12455, "Random number generator seed");
      29          23 :   return params;
      30           0 : }
      31             : 
      32          12 : FourierNoise::FourierNoise(const InputParameters & parameters)
      33             :   : Function(parameters),
      34          12 :     _lambda(getParam<Real>("lambda")),
      35          36 :     _fe_problem(*getCheckedPointerParam<FEProblemBase *>("_fe_problem_base"))
      36             : {
      37             :   MooseRandom rng;
      38          24 :   rng.seed(0, getParam<unsigned int>("seed"));
      39             : 
      40          24 :   if (isParamValid("num_terms"))
      41             :   {
      42             :     // random terms
      43           0 :     _series.resize(getParam<unsigned int>("num_terms"));
      44           0 :     const Real scale = 2.0 * (2.0 * libMesh::pi) / _lambda;
      45             : 
      46             :     // check
      47           0 :     if (_series.empty())
      48           0 :       paramError("num_terms",
      49             :                  "If specifying the number of terms, supply a number greater than zero.");
      50             : 
      51             :     // fill terms
      52           0 :     for (auto & f : _series)
      53             :     {
      54             :       // get a vector with length <= 0.5
      55             :       Real r2;
      56             :       do
      57             :       {
      58           0 :         const Real x = rng.rand(0) - 0.5;
      59           0 :         const Real y = rng.rand(0) - 0.5;
      60           0 :         f.k = RealVectorValue(x, y, 0.0);
      61             :         r2 = f.k.norm_sq();
      62           0 :       } while (r2 > 0.25);
      63             : 
      64             :       // scale maximum to a wavelength of lambda
      65             :       f.k *= scale;
      66             : 
      67           0 :       f.c = rng.randNormal(0, 0.0, 1.0);
      68           0 :       f.s = rng.randNormal(0, 0.0, 1.0);
      69             :     }
      70             :   }
      71             :   else
      72             :   {
      73             :     // k-space grid resulting in periodic noise
      74          12 :     MooseMesh & mesh = _fe_problem.mesh();
      75          12 :     if (!mesh.isRegularOrthogonal())
      76           0 :       mooseError("Periodic Fourier Noise requires a regular orthogonal mesh.");
      77             : 
      78          12 :     const Real dx = 2.0 * libMesh::pi / mesh.dimensionWidth(0);
      79          12 :     const Real dy = 2.0 * libMesh::pi / mesh.dimensionWidth(1);
      80          12 :     const Real rmax = 2.0 * libMesh::pi / _lambda;
      81             : 
      82          12 :     const int xmax = rmax / dx;
      83          12 :     const int ymax = rmax / dy;
      84             : 
      85          12 :     const Real rmax2 = rmax * rmax;
      86          72 :     for (int x = 0; x < xmax; ++x)
      87         660 :       for (int y = -ymax; y < ymax; ++y)
      88         600 :         if (x > 0 || y > 0)
      89             :         {
      90             :           SeriesItem f;
      91         528 :           f.k = RealVectorValue(x * dx, y * dy, 0.0);
      92         528 :           if (f.k.norm_sq() <= rmax2)
      93             :           {
      94         456 :             f.c = rng.randNormal(0, 0.0, 1.0);
      95         456 :             f.s = rng.randNormal(0, 0.0, 1.0);
      96         456 :             _series.push_back(f);
      97             :           }
      98             :         }
      99             :   }
     100             : 
     101          12 :   _scale = std::sqrt(1.0 / _series.size());
     102          12 : }
     103             : 
     104             : Real
     105        9600 : FourierNoise::value(Real, const Point & p) const
     106             : {
     107             :   Real v = 0.0;
     108      374400 :   for (const auto & f : _series)
     109      364800 :     v += f.s * std::sin(p * f.k) + f.c * std::cos(p * f.k);
     110        9600 :   return v * _scale;
     111             : }
     112             : 
     113             : ADReal
     114           0 : FourierNoise::value(const ADReal &, const ADPoint & p) const
     115             : {
     116           0 :   ADReal v = 0.0;
     117           0 :   for (const auto & f : _series)
     118           0 :     v += f.s * std::sin(p * f.k) + f.c * std::cos(p * f.k);
     119           0 :   return v * _scale;
     120             : }