Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
       3             : 
       4             : // This library is free software; you can redistribute it and/or
       5             : // modify it under the terms of the GNU Lesser General Public
       6             : // License as published by the Free Software Foundation; either
       7             : // version 2.1 of the License, or (at your option) any later version.
       8             : 
       9             : // This library is distributed in the hope that it will be useful,
      10             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      11             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12             : // Lesser General Public License for more details.
      13             : 
      14             : // You should have received a copy of the GNU Lesser General Public
      15             : // License along with this library; if not, write to the Free Software
      16             : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
      17             : 
      18             : #include "libmesh/adaptive_time_solver.h"
      19             : #include "libmesh/diff_system.h"
      20             : #include "libmesh/dof_map.h"
      21             : #include "libmesh/numeric_vector.h"
      22             : #include "libmesh/no_solution_history.h"
      23             : 
      24             : // Debug
      25             : #include "libmesh/system_norm.h"
      26             : #include "libmesh/enum_norm_type.h"
      27             : 
      28             : namespace libMesh
      29             : {
      30             : 
      31             : 
      32             : 
      33         420 : AdaptiveTimeSolver::AdaptiveTimeSolver (sys_type & s)
      34             :   : FirstOrderUnsteadySolver(s),
      35         396 :     core_time_solver(),
      36         396 :     target_tolerance(1.e-3),
      37         396 :     upper_tolerance(0.0),
      38         396 :     max_deltat(0.),
      39         396 :     min_deltat(0.),
      40         396 :     max_growth(0.),
      41         420 :     completed_timestep_size(s.deltat),
      42         432 :     global_tolerance(true)
      43             : {
      44             :   // the child class must populate core_time_solver
      45             :   // with whatever actual time solver is to be used
      46         420 : }
      47             : 
      48             : 
      49             : 
      50         396 : AdaptiveTimeSolver::~AdaptiveTimeSolver () = default;
      51             : 
      52             : 
      53             : 
      54         420 : void AdaptiveTimeSolver::init()
      55             : {
      56          12 :   libmesh_assert(core_time_solver.get());
      57             : 
      58             :   // We override this because our core_time_solver is the one that
      59             :   // needs to handle new vectors, diff_solver->init(), etc
      60         420 :   core_time_solver->init();
      61             : 
      62             :   // Set the core_time_solver's solution history object to be the same one as
      63             :   // that for the outer adaptive time solver
      64         420 :   core_time_solver->set_solution_history((this->get_solution_history()));
      65             : 
      66             :   // Now that we have set the SolutionHistory object for the coretimesolver,
      67             :   // we set the SolutionHistory type for the timesolver to be NoSolutionHistory
      68             :   // All storage and retrieval will be handled by the coretimesolver directly.
      69          24 :   NoSolutionHistory outersolver_solution_history;
      70         420 :   this->set_solution_history(outersolver_solution_history);
      71             : 
      72             :   // As an UnsteadySolver, we have an old_local_nonlinear_solution, but it
      73             :   // isn't pointing to the right place - fix it
      74          12 :   old_local_nonlinear_solution = core_time_solver->old_local_nonlinear_solution;
      75         420 : }
      76             : 
      77             : 
      78             : 
      79         560 : void AdaptiveTimeSolver::reinit()
      80             : {
      81          16 :   libmesh_assert(core_time_solver.get());
      82             : 
      83             :   // We override this because our core_time_solver is the one that
      84             :   // needs to handle new vectors, diff_solver->reinit(), etc
      85         560 :   core_time_solver->reinit();
      86         560 : }
      87             : 
      88             : 
      89        4200 : void AdaptiveTimeSolver::advance_timestep ()
      90             : {
      91             :   // The first access of advance_timestep happens via solve, not user code
      92             :   // It is used here to store any initial conditions data
      93        4200 :   if (!first_solve)
      94             :     {
      95        3780 :       _system.time += this->completed_timestep_size;
      96             :     }
      97             :   else
      98             :     {
      99             :     // We are here because of a call to advance_timestep that happens
     100             :     // via solve, the very first solve. All we are doing here is storing
     101             :     // the initial condition. The actual solution computed via this solve
     102             :     // will be stored when we call advance_timestep in the user's timestep loop
     103         420 :     first_solve = false;
     104          12 :     core_time_solver->set_first_solve(false);
     105             :     }
     106             : 
     107             :   // For the adaptive time solver, all SH operations
     108             :   // are handled by the core_time_solver's SH object
     109             :   // Sub solution storage is handled internally by the core time solver,
     110             :   // but the 'full step' solution is stored here to maintain consistency
     111             :   // with the fixed timestep scheme.
     112        4200 :   core_time_solver->get_solution_history().store(false, _system.time);
     113             : 
     114             :   NumericVector<Number> & old_nonlinear_soln =
     115        4200 :     _system.get_vector("_old_nonlinear_solution");
     116             :   NumericVector<Number> & nonlinear_solution =
     117        4200 :     *(_system.solution);
     118             : 
     119        4200 :   old_nonlinear_soln = nonlinear_solution;
     120             : 
     121             :   old_nonlinear_soln.localize
     122        4200 :     (*old_local_nonlinear_solution,
     123        4200 :      _system.get_dof_map().get_send_list());
     124        4200 : }
     125             : 
     126        4200 : void AdaptiveTimeSolver::adjoint_advance_timestep ()
     127             : {
     128             :   // Store the computed full step adjoint solution for future use (sub steps are handled internally by the core time solver)
     129        4200 :   core_time_solver->get_solution_history().store(true, _system.time);
     130             : 
     131             :   // For the first adjoint step ensure that we use the last primal timestep.
     132        4200 :   if (first_adjoint_step)
     133             :   {
     134         432 :     _system.deltat = dynamic_cast<DifferentiableSystem &>(_system).time_solver->TimeSolver::last_completed_timestep_size();
     135         420 :     first_adjoint_step = false;
     136             :   }
     137             : 
     138             :   // Before moving to the next time instant, copy over the current adjoint solutions into _old_adjoint_solutions
     139       11480 :   for(auto i : make_range(_system.n_qois()))
     140             :   {
     141        7488 :    std::string old_adjoint_solution_name = "_old_adjoint_solution";
     142        7072 :    old_adjoint_solution_name+= std::to_string(i);
     143        7488 :    NumericVector<Number> & old_adjoint_solution_i = _system.get_vector(old_adjoint_solution_name);
     144        7280 :    NumericVector<Number> & adjoint_solution_i = _system.get_adjoint_solution(i);
     145        7280 :    old_adjoint_solution_i = adjoint_solution_i;
     146             :   }
     147             : 
     148        4200 :   _system.time -= _system.deltat;
     149             : 
     150             :   // For the adaptive time solver, all SH operations
     151             :   // are handled by the core_time_solver's SH object
     152             :   // Retrieve the primal solution for the next adjoint calculation,
     153             :   // by using the core time solver's solution history object.
     154        4200 :   core_time_solver->get_solution_history().retrieve(true, _system.time);
     155             : 
     156             :   // We also need to tell the core time solver that the adjoint initial conditions have been set
     157         120 :   core_time_solver->set_first_adjoint_step(false);
     158             : 
     159             :   // Dont forget to localize the old_nonlinear_solution !
     160        4200 :   _system.get_vector("_old_nonlinear_solution").localize
     161        4200 :     (*old_local_nonlinear_solution,
     162        4200 :     _system.get_dof_map().get_send_list());
     163        4200 : }
     164             : 
     165         700 : void AdaptiveTimeSolver::retrieve_timestep ()
     166             : {
     167             :   // Ask the core time solver to retrieve all the stored vectors
     168             :   // at the current time
     169         700 :   core_time_solver->retrieve_timestep();
     170         700 : }
     171             : 
     172           0 : Real AdaptiveTimeSolver::error_order () const
     173             : {
     174           0 :   libmesh_assert(core_time_solver.get());
     175             : 
     176           0 :   return core_time_solver->error_order();
     177             : }
     178             : 
     179             : 
     180             : 
     181    10128536 : bool AdaptiveTimeSolver::element_residual (bool request_jacobian,
     182             :                                            DiffContext & context)
     183             : {
     184      920776 :   libmesh_assert(core_time_solver.get());
     185             : 
     186    10128536 :   return core_time_solver->element_residual(request_jacobian, context);
     187             : }
     188             : 
     189             : 
     190             : 
     191     1277760 : bool AdaptiveTimeSolver::side_residual (bool request_jacobian,
     192             :                                         DiffContext & context)
     193             : {
     194       63888 :   libmesh_assert(core_time_solver.get());
     195             : 
     196      702768 :   return core_time_solver->side_residual(request_jacobian, context);
     197             : }
     198             : 
     199             : 
     200             : 
     201           0 : bool AdaptiveTimeSolver::nonlocal_residual (bool request_jacobian,
     202             :                                             DiffContext & context)
     203             : {
     204           0 :   libmesh_assert(core_time_solver.get());
     205             : 
     206           0 :   return core_time_solver->nonlocal_residual(request_jacobian, context);
     207             : }
     208             : 
     209             : 
     210             : 
     211       16660 : std::unique_ptr<DiffSolver> & AdaptiveTimeSolver::diff_solver()
     212             : {
     213       16660 :   return core_time_solver->diff_solver();
     214             : }
     215             : 
     216             : 
     217             : 
     218      101360 : std::unique_ptr<LinearSolver<Number>> & AdaptiveTimeSolver::linear_solver()
     219             : {
     220      101360 :   return core_time_solver->linear_solver();
     221             : }
     222             : 
     223             : 
     224             : 
     225       11340 : Real AdaptiveTimeSolver::calculate_norm(System & s,
     226             :                                         NumericVector<Number> & v)
     227             : {
     228       11340 :   return s.calculate_norm(v, component_norm);
     229             : }
     230             : 
     231             : } // namespace libMesh