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 "AB2PredictorCorrector.h"
      11             : #include "AdamsPredictor.h"
      12             : #include "Problem.h"
      13             : #include "FEProblem.h"
      14             : #include "MooseApp.h"
      15             : #include "NonlinearSystem.h"
      16             : #include "AuxiliarySystem.h"
      17             : #include "TimeIntegrator.h"
      18             : #include "Conversion.h"
      19             : 
      20             : #include "libmesh/nonlinear_solver.h"
      21             : #include "libmesh/numeric_vector.h"
      22             : 
      23             : // C++ Includes
      24             : #include <iomanip>
      25             : #include <iostream>
      26             : #include <fstream>
      27             : 
      28             : using namespace libMesh;
      29             : 
      30             : registerMooseObject("MooseApp", AB2PredictorCorrector);
      31             : 
      32             : InputParameters
      33        3105 : AB2PredictorCorrector::validParams()
      34             : {
      35        3105 :   InputParameters params = TimeStepper::validParams();
      36        6210 :   params.addClassDescription(
      37             :       "Implements second order Adams-Bashforth method for timestep calculation.");
      38       12420 :   params.addRequiredParam<Real>("e_tol", "Target error tolerance.");
      39       12420 :   params.addRequiredParam<Real>("e_max", "Maximum acceptable error.");
      40       12420 :   params.addRequiredParam<Real>("dt", "Initial time step size");
      41       12420 :   params.addParam<Real>("max_increase", 1.0e9, "Maximum ratio that the time step can increase.");
      42        9315 :   params.addParam<int>(
      43        6210 :       "steps_between_increase", 1, "the number of time steps before recalculating dt");
      44       12420 :   params.addParam<int>("start_adapting", 2, "when to start taking adaptive time steps");
      45        9315 :   params.addParam<Real>("scaling_parameter", .8, "scaling parameter for dt selection");
      46        3105 :   return params;
      47           0 : }
      48             : 
      49          22 : AB2PredictorCorrector::AB2PredictorCorrector(const InputParameters & parameters)
      50             :   : TimeStepper(parameters),
      51          44 :     _u1(_fe_problem.getNonlinearSystemBase(/*nl_sys=*/0).addVector("u1", true, GHOSTED)),
      52          44 :     _aux1(_fe_problem.getAuxiliarySystem().addVector("aux1", true, GHOSTED)),
      53          44 :     _pred1(_fe_problem.getNonlinearSystemBase(/*nl_sys=*/0).addVector("pred1", true, GHOSTED)),
      54          44 :     _dt_full(declareRestartableData<Real>("dt_full", 0)),
      55          44 :     _error(declareRestartableData<Real>("error", 0)),
      56          44 :     _e_tol(getParam<Real>("e_tol")),
      57          44 :     _e_max(getParam<Real>("e_max")),
      58          44 :     _max_increase(getParam<Real>("max_increase")),
      59          44 :     _steps_between_increase(getParam<int>("steps_between_increase")),
      60          44 :     _dt_steps_taken(declareRestartableData<int>("dt_steps_taken", 0)),
      61          44 :     _start_adapting(getParam<int>("start_adapting")),
      62          44 :     _my_dt_old(declareRestartableData<Real>("my_dt_old", 0)),
      63          44 :     _infnorm(declareRestartableData<Real>("infnorm", 0)),
      64          88 :     _scaling_parameter(getParam<Real>("scaling_parameter"))
      65             : {
      66          22 :   Real predscale = 1.;
      67          44 :   InputParameters params = _app.getFactory().getValidParams("AdamsPredictor");
      68          44 :   params.set<Real>("scale") = predscale;
      69          88 :   _fe_problem.addPredictor("AdamsPredictor", "adamspredictor", params);
      70          22 : }
      71             : 
      72             : void
      73          22 : AB2PredictorCorrector::preExecute()
      74             : {
      75          22 :   TimeStepper::preExecute();
      76          22 : }
      77             : 
      78             : void
      79          80 : AB2PredictorCorrector::preSolve()
      80             : {
      81             :   // save dt
      82          80 :   _dt_full = _dt;
      83          80 : }
      84             : 
      85             : void
      86          80 : AB2PredictorCorrector::step()
      87             : {
      88          80 :   NonlinearSystemBase & nl = _fe_problem.getNonlinearSystemBase(/*nl_sys=*/0);
      89          80 :   AuxiliarySystem & aux = _fe_problem.getAuxiliarySystem();
      90             : 
      91          80 :   _fe_problem.solve(/*nl_sys=*/0);
      92          80 :   _converged = _fe_problem.converged(/*nl_sys=*/0);
      93          80 :   if (_converged)
      94             :   {
      95          80 :     _u1 = *nl.currentSolution();
      96          80 :     _u1.close();
      97             : 
      98          80 :     _aux1 = *aux.currentSolution();
      99          80 :     _aux1.close();
     100          80 :     if (_t_step >= _start_adapting)
     101             :     {
     102             :       // Calculate error if past the first solve
     103          60 :       _error = estimateTimeError(_u1);
     104             : 
     105          60 :       _infnorm = _u1.linfty_norm();
     106          60 :       _e_max = 1.1 * _e_tol * _infnorm;
     107          60 :       _console << "Time Error Estimate: " << _error << std::endl;
     108             :     }
     109             :     else
     110             :     {
     111             :       // First time step is problematic, sure we converged but what does that mean? We don't know.
     112             :       // Nor can we calculate the error on the first time step.
     113             :     }
     114             :   }
     115          80 : }
     116             : 
     117             : bool
     118         339 : AB2PredictorCorrector::converged() const
     119             : {
     120         339 :   if (!_converged)
     121           0 :     return false;
     122         339 :   if (_error < _e_max)
     123         339 :     return true;
     124             :   else
     125           0 :     return false;
     126             : }
     127             : 
     128             : void
     129          80 : AB2PredictorCorrector::postSolve()
     130             : {
     131          80 :   if (!converged())
     132           0 :     _dt_steps_taken = 0;
     133             : 
     134          80 :   if (_error >= _e_max)
     135           0 :     _console << "Marking last solve not converged " << _error << " " << _e_max << std::endl;
     136          80 : }
     137             : 
     138             : Real
     139          60 : AB2PredictorCorrector::computeDT()
     140             : {
     141          60 :   if (_t_step <= _start_adapting)
     142          20 :     return _dt;
     143             : 
     144          40 :   _my_dt_old = _dt;
     145             : 
     146          40 :   _dt_steps_taken += 1;
     147          40 :   if (_dt_steps_taken >= _steps_between_increase)
     148             :   {
     149             : 
     150          40 :     Real new_dt = _dt_full * _scaling_parameter * std::pow(_infnorm * _e_tol / _error, 1.0 / 3.0);
     151             : 
     152          40 :     if (new_dt / _dt_full > _max_increase)
     153           0 :       new_dt = _dt_full * _max_increase;
     154          40 :     _dt_steps_taken = 0;
     155          40 :     return new_dt;
     156             :   }
     157             : 
     158           0 :   return _dt;
     159             : }
     160             : 
     161             : Real
     162          40 : AB2PredictorCorrector::computeInitialDT()
     163             : {
     164         120 :   return getParam<Real>("dt");
     165             : }
     166             : 
     167             : Real
     168          60 : AB2PredictorCorrector::estimateTimeError(NumericVector<Number> & solution)
     169             : {
     170          60 :   _pred1 = _fe_problem.getNonlinearSystemBase(/*nl_sys=*/0).getPredictor()->solutionPredictor();
     171             :   const auto & ti =
     172          60 :       _fe_problem.getNonlinearSystemBase(/*nl_sys=*/0).getTimeIntegrator(/*var_num=*/0);
     173             : 
     174          60 :   auto scheme = Moose::stringToEnum<Moose::TimeIntegratorType>(ti.type());
     175          60 :   Real dt_old = _my_dt_old;
     176          60 :   if (dt_old == 0)
     177          20 :     dt_old = _dt;
     178             : 
     179          60 :   switch (scheme)
     180             :   {
     181           0 :     case Moose::TI_IMPLICIT_EULER:
     182             :     {
     183           0 :       _pred1 *= -1;
     184           0 :       _pred1 += solution;
     185           0 :       Real calc = _dt * _dt * .5;
     186           0 :       _pred1 *= calc;
     187           0 :       return _pred1.l2_norm();
     188             :     }
     189           0 :     case Moose::TI_CRANK_NICOLSON:
     190             :     {
     191           0 :       _pred1 -= solution;
     192           0 :       _pred1 *= (_dt) / (3.0 * (_dt + dt_old));
     193           0 :       return _pred1.l2_norm();
     194             :     }
     195          60 :     case Moose::TI_BDF2:
     196             :     {
     197          60 :       _pred1 *= -1.0;
     198          60 :       _pred1 += solution;
     199          60 :       Real topcalc = 2.0 * (_dt + dt_old) * (_dt + dt_old);
     200          60 :       Real bottomcalc = 6.0 * _dt * _dt + 12.0 * _dt * dt_old + 5.0 * dt_old * dt_old;
     201          60 :       _pred1 *= topcalc / bottomcalc;
     202             : 
     203          60 :       return _pred1.l2_norm();
     204             :     }
     205           0 :     default:
     206           0 :       break;
     207             :   }
     208           0 :   return -1;
     209             : }