Line data    Source code

       1             : /**********************************************************************/
       2             : /*                     DO NOT MODIFY THIS HEADER                      */
       3             : /* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
       4             : /*                                                                    */
       5             : /*            Copyright 2017 Battelle Energy Alliance, LLC            */
       6             : /*                        ALL RIGHTS RESERVED                         */
       7             : /**********************************************************************/
       8             : 
       9             : #include "SpectralExecutionerDiffusion.h"
      10             : #include "InputParameters.h"
      11             : 
      12             : registerMooseObject("MagpieApp", SpectralExecutionerDiffusion);
      13             : 
      14             : InputParameters
      15           6 : SpectralExecutionerDiffusion::validParams()
      16             : {
      17           6 :   InputParameters params = SpectralExecutionerBase::validParams();
      18           6 :   params.addClassDescription("Executioner for spectral solve of diffusion equation");
      19          12 :   params.addParam<Real>("diffusion_coefficient", 1.0, "Diffusion coefficient");
      20          12 :   params.addParam<Real>("time_step", 1.0, "Time step for ODE integration");
      21          12 :   params.addParam<unsigned int>("number_steps", 0.0, "Time step for ODE integration");
      22           6 :   return params;
      23           0 : }
      24             : 
      25           3 : SpectralExecutionerDiffusion::SpectralExecutionerDiffusion(const InputParameters & parameters)
      26             :   : SpectralExecutionerBase(parameters),
      27           3 :     _diff_coeff(getParam<Real>("diffusion_coefficient")),
      28           6 :     _dt(getParam<Real>("time_step")),
      29           9 :     _nsteps(getParam<unsigned int>("number_steps"))
      30             : {
      31           3 :   _t_current = 0.0;
      32           3 : }
      33             : 
      34             : void
      35           0 : SpectralExecutionerDiffusion::executeExplicit()
      36             : {
      37             :   unsigned int thisStep = 0;
      38             : 
      39           0 :   _time_step = thisStep;
      40           0 :   _time = _time_step;
      41           0 :   _fe_problem.outputStep(EXEC_INITIAL);
      42           0 :   _fe_problem.advanceState();
      43             : 
      44           0 :   auto & u_buffer = getFFTBuffer<Real>("u");
      45             :   auto u = u_buffer.realSpace();
      46             : 
      47           0 :   for (unsigned int step_no = 0; step_no < _nsteps; step_no++)
      48             :   {
      49             :     // Forward transform to get \hat{u}_{k}
      50           0 :     u_buffer.forwardRaw();
      51             : 
      52           0 :     advanceDiffusionStep(u_buffer, _diff_coeff);
      53             : 
      54           0 :     u_buffer.backward();
      55             : 
      56           0 :     thisStep++;
      57           0 :     _t_current += _dt;
      58           0 :     _time_step = thisStep;
      59             : 
      60           0 :     _fe_problem.execute(EXEC_FINAL);
      61           0 :     _time = _t_current;
      62             : 
      63             :     Moose::out << "_t_current: " << _t_current << ". \n";
      64             : 
      65           0 :     _fe_problem.outputStep(EXEC_FINAL);
      66             : 
      67           0 :     if (step_no != _nsteps - 1)
      68           0 :       _fe_problem.advanceState();
      69             :   }
      70           0 : }
      71             : 
      72             : void
      73           3 : SpectralExecutionerDiffusion::getGreensFunction(FFTBufferBase<Real> & greens,
      74             :                                                 Real time,
      75             :                                                 const Real D)
      76             : {
      77           3 :   Real accGreens = 0.0;
      78             : 
      79             :   const Point & box_size = greens.getBoxSize();
      80             : 
      81             :   const auto & grid = greens.grid();
      82           3 :   switch (greens.dim())
      83             :   {
      84           0 :     case 1:
      85             :     {
      86           0 :       mooseError("Error: Dimension 1 not implemented yet.");
      87             :       break;
      88             :     }
      89             : 
      90           3 :     case 2:
      91             :     {
      92             :       std::size_t index = 0;
      93             : 
      94           3 :       const int ni = grid[0];
      95           3 :       const int nj = grid[1]; // Real space.
      96             : 
      97          63 :       for (int i = 0; i < ni; ++i)
      98        1260 :         for (int j = 0; j < nj; ++j)
      99             :         {
     100             : 
     101        1200 :           Real x1 = Real(i) / Real(ni) * box_size(0);
     102        1200 :           Real y1 = Real(j) / Real(nj) * box_size(1);
     103             : 
     104        1200 :           Real x2 = box_size(0) - x1;
     105        1200 :           Real y2 = box_size(1) - y1;
     106             : 
     107             :           Moose::out << "box size: " << box_size(0) << ", " << box_size(1) << "\n";
     108             :           Moose::out << "number of cells: " << ni << ", " << nj << "\n";
     109             : 
     110        1200 :           if (time == 0)
     111           0 :             mooseError("Greens function undefined at t = 0s in the FFT solver.");
     112             :           else
     113             :           {
     114        1200 :             Real sum = std::exp(-(x1 * x1 + y1 * y1) / 4 / D / time) +
     115        1200 :                        std::exp(-(x1 * x1 + y2 * y2) / 4 / D / time) +
     116        1200 :                        std::exp(-(x2 * x2 + y1 * y1) / 4 / D / time) +
     117        1200 :                        std::exp(-(x2 * x2 + y2 * y2) / 4 / D / time);
     118             : 
     119        1200 :             if (sum == 0 && i != 0 && j != 0)
     120           0 :               mooseError("Precision lost in the analytical integration of heat equation");
     121             : 
     122        1200 :             Real correction_factor = (Real(ni) / box_size(0)) * (Real(nj) / box_size(1));
     123             : 
     124        1200 :             greens.realSpace()[index] =
     125        1200 :                 (1.0 / (4 * libMesh::pi * D * time)) * sum / correction_factor;
     126             :           }
     127        1200 :           accGreens += greens.realSpace()[index];
     128        1200 :           index++;
     129             :         }
     130           3 :       mooseInfo("Accumulated value of Greens function: ", accGreens, "\n");
     131             :       break;
     132             :     }
     133             : 
     134           0 :     case 3:
     135             :     {
     136           0 :       mooseError("Error: Dimension 3 not implemented yet.");
     137             :       break;
     138             :     }
     139             : 
     140           0 :     default:
     141           0 :       mooseError("Invalid buffer dimension");
     142             :   }
     143           3 : }
     144             : void
     145           3 : SpectralExecutionerDiffusion::execute()
     146             : {
     147             :   unsigned int thisStep = 0;
     148             : 
     149           3 :   _time_step = thisStep;
     150           3 :   _time = _time_step;
     151           3 :   _fe_problem.outputStep(EXEC_INITIAL);
     152           3 :   _fe_problem.advanceState();
     153             : 
     154           3 :   auto & u_buffer = getFFTBuffer<Real>("u");
     155           3 :   auto & greens = getFFTBuffer<Real>("greens_buffer");
     156             : 
     157             :   // Get Greens function of
     158           3 :   getGreensFunction(greens, _dt, _diff_coeff);
     159           3 :   greens.forwardRaw();
     160             : 
     161          63 :   for (unsigned int step_no = 0; step_no < _nsteps; step_no++)
     162             :   {
     163          60 :     u_buffer.forwardRaw();
     164             : 
     165          60 :     u_buffer.reciprocalSpace() *= greens.reciprocalSpace();
     166             :     // scalarMultiplyBuffer(u_buffer, greens);
     167          60 :     u_buffer.backward();
     168             : 
     169             :     // End of diffusion computations
     170          60 :     thisStep++;
     171          60 :     _t_current += _dt;
     172          60 :     _time_step = thisStep;
     173             : 
     174          60 :     _fe_problem.execute(EXEC_FINAL);
     175          60 :     _time = _t_current;
     176             : 
     177             :     Moose::out << "_t_current: " << _t_current << ". \n";
     178             : 
     179          60 :     _fe_problem.outputStep(EXEC_FINAL);
     180             : 
     181          60 :     if (step_no != _nsteps - 1)
     182          57 :       _fe_problem.advanceState();
     183             :   }
     184           3 : }
     185             : 
     186             : void
     187           0 : SpectralExecutionerDiffusion::executeTotalDiffusion()
     188             : {
     189             :   unsigned int thisStep = 0;
     190             : 
     191           0 :   _time_step = thisStep;
     192           0 :   _time = _time_step;
     193           0 :   _fe_problem.outputStep(EXEC_INITIAL);
     194           0 :   _fe_problem.advanceState();
     195             : 
     196           0 :   FFTBufferBase<Real> & u_buffer = getFFTBuffer<Real>("u");
     197           0 :   FFTBufferBase<Real> & u_init_buffer = getFFTBuffer<Real>("u_initial");
     198             : 
     199           0 :   u_init_buffer.forwardRaw();
     200             :   // auto &greens = getFFTBuffer<Real>("greens_buffer");
     201             : 
     202           0 :   for (unsigned int step_no = 0; step_no < _nsteps; step_no++)
     203             :   {
     204             :     // Forward transform to get \hat{u}_{k}
     205           0 :     u_buffer.forwardRaw();
     206             : 
     207           0 :     advanceDiffusionStepTotal(u_init_buffer, u_buffer, _diff_coeff, _time);
     208             : 
     209           0 :     u_buffer.backward();
     210             : 
     211           0 :     thisStep++;
     212           0 :     _t_current += _dt;
     213           0 :     _time_step = thisStep;
     214             : 
     215           0 :     _fe_problem.execute(EXEC_FINAL);
     216           0 :     _time = _t_current;
     217             : 
     218             :     Moose::out << "_t_current: " << _t_current << ". \n";
     219             : 
     220           0 :     _fe_problem.outputStep(EXEC_FINAL);
     221             : 
     222           0 :     if (step_no != _nsteps - 1)
     223           0 :       _fe_problem.advanceState();
     224             :   }
     225           0 : }
     226             : 
     227             : void
     228           0 : SpectralExecutionerDiffusion::scalarMultiplyBuffer(FFTBufferBase<Real> & u,
     229             :                                                    const FFTBufferBase<Real> & greens)
     230             : {
     231             :   const auto & grid = u.grid();
     232             : 
     233           0 :   switch (u.dim())
     234             :   {
     235           0 :     case 1:
     236             :     {
     237           0 :       mooseError("Error: Dimension 1 not implemented yet.");
     238             :       return;
     239             :     }
     240             : 
     241           0 :     case 2:
     242             :     {
     243             :       std::size_t index = 0;
     244           0 :       const int ni = grid[0];
     245           0 :       const int nj = (grid[1] >> 1) + 1;
     246             : 
     247           0 :       for (int i = 0; i < ni; ++i)
     248           0 :         for (int j = 0; j < nj; ++j)
     249             :         {
     250           0 :           u.reciprocalSpace()[index] = greens.reciprocalSpace()[index] * u.reciprocalSpace()[index];
     251           0 :           index++;
     252             :         }
     253           0 :       return;
     254             :     }
     255             : 
     256           0 :     case 3:
     257             :     {
     258           0 :       mooseError("Error: Dimension 3 not implemented yet.");
     259             :       return;
     260             :     }
     261             : 
     262           0 :     default:
     263           0 :       mooseError("Invalid buffer dimension");
     264             :   }
     265             : }
     266             : 
     267             : void
     268           0 : SpectralExecutionerDiffusion::advanceDiffusionStepTotal(const FFTBufferBase<Real> & u_initial,
     269             :                                                         FFTBufferBase<Real> & u,
     270             :                                                         const Real D,
     271             :                                                         const Real time)
     272             : {
     273             :   const auto & grid = u.grid();
     274           0 :   switch (u.dim())
     275             :   {
     276           0 :     case 1:
     277             :     {
     278           0 :       mooseError("Error: Dimension 1 not implemented yet.");
     279             :       return;
     280             :     }
     281             : 
     282           0 :     case 2:
     283             :     {
     284             :       std::size_t index = 0;
     285           0 :       const int ni = grid[0];
     286           0 :       const int nj = (grid[1] >> 1) + 1;
     287             : 
     288             :       const auto & ivec = u.kTable(0);
     289             :       const auto & jvec = u.kTable(1);
     290             : 
     291           0 :       for (int i = 0; i < ni; ++i)
     292           0 :         for (int j = 0; j < nj; ++j)
     293             :         {
     294           0 :           u.reciprocalSpace()[index] =
     295             :               u_initial.reciprocalSpace()[index] *
     296           0 :               std::exp(-D * (ivec[i] * ivec[i] + jvec[j] * jvec[j]) * time);
     297           0 :           index++;
     298             :         }
     299           0 :       return;
     300             :     }
     301             : 
     302           0 :     case 3:
     303             :     {
     304           0 :       mooseError("Error: Dimension 3 not implemented yet.");
     305             :       return;
     306             :     }
     307             : 
     308           0 :     default:
     309           0 :       mooseError("Invalid buffer dimension");
     310             :   }
     311             : }
     312             : 
     313             : void
     314           0 : SpectralExecutionerDiffusion::advanceDiffusionStep(FFTBufferBase<Real> & inOut, const Real D)
     315             : {
     316             :   const Real sizePixel = 1.0; // Assumed same size in all directions
     317             : 
     318             :   const auto & grid = inOut.grid();
     319           0 :   switch (inOut.dim())
     320             :   {
     321           0 :     case 1:
     322             :     {
     323           0 :       mooseError("Error: Dimension 1 not implemented yet.");
     324             :       return;
     325             :     }
     326             : 
     327           0 :     case 2:
     328             :     {
     329             :       std::size_t index = 0;
     330           0 :       const int ni = grid[0];
     331           0 :       const int nj = (grid[1] >> 1) + 1;
     332             : 
     333           0 :       for (int i = 0; i < ni; ++i)
     334           0 :         for (int j = 0; j < nj; ++j)
     335             :         {
     336           0 :           Real r1 = D * _dt / (sizePixel * sizePixel);
     337             :           Real r2 = D * _dt / (sizePixel * sizePixel);
     338             : 
     339           0 :           Real hk1 = (1 - 2 * r1 * (1 - std::cos(libMesh::pi * i / ((ni) / 2))));
     340           0 :           Real hk2 = (1 - 2 * r2 * (1 - std::cos(libMesh::pi * j / ((nj) / 2))));
     341           0 :           inOut.reciprocalSpace()[index] *= hk1 * hk2;
     342           0 :           index++;
     343             :         }
     344           0 :       return;
     345             :     }
     346             : 
     347           0 :     case 3:
     348             :     {
     349           0 :       mooseError("Error: Dimension 3 not implemented yet.");
     350             :       return;
     351             :     }
     352             : 
     353           0 :     default:
     354           0 :       mooseError("Invalid buffer dimension");
     355             :   }
     356             : }