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             : // Local includes
      19             : #include "libmesh/memory_solution_history.h"
      20             : #include "libmesh/memory_history_data.h"
      21             : 
      22             : #include "libmesh/diff_system.h"
      23             : 
      24             : #include <cmath>
      25             : #include <iterator>
      26             : 
      27             : namespace libMesh
      28             : {
      29             : 
      30        2240 : MemorySolutionHistory::~MemorySolutionHistory ()
      31             : {
      32        2240 : }
      33             : 
      34             : // This functions saves all the projected system vectors for
      35             : // future use
      36       17080 : void MemorySolutionHistory::store(bool /* is_adjoint_solve */, Real time)
      37             : {
      38       17080 :   this->find_stored_entry(time, true);
      39             : 
      40             :   // In an empty history we create the first entry
      41       17080 :   if (stored_data.begin() == stored_data.end())
      42             :     {
      43         576 :       stored_data[time] = std::make_unique<MemoryHistoryData>(_system);
      44         560 :       stored_datum = stored_data.begin();
      45             :     }
      46             : 
      47             :   // If we're past the end we can create a new entry
      48       17080 :   if (time - stored_datum->first > TOLERANCE )
      49             :     {
      50             : #ifndef NDEBUG
      51         228 :       ++stored_datum;
      52         228 :       libmesh_assert (stored_datum == stored_data.end());
      53             : #endif
      54        7980 :       stored_data[time] = std::make_unique<MemoryHistoryData>(_system);
      55        7980 :       stored_datum = stored_data.end();
      56         228 :       --stored_datum;
      57             :     }
      58             : 
      59             :   // If we're before the beginning we can create a new entry
      60        9100 :   else if (stored_datum->first - time > TOLERANCE)
      61             :     {
      62           0 :       libmesh_assert (stored_datum == stored_data.begin());
      63           0 :       stored_data[time] = std::make_unique<MemoryHistoryData>(_system);
      64           0 :       stored_datum = stored_data.begin();
      65             :     }
      66             : 
      67             :   // We don't support inserting entries elsewhere
      68         488 :   libmesh_assert(std::abs(stored_datum->first - time) < TOLERANCE);
      69             : 
      70             :   // First we handle the case of the initial data, this is the only case in which
      71             :   // stored_data will have size one
      72       17080 :   if(stored_data.size() == 1)
      73             :   {
      74             :     // The initial data should only be stored once.
      75         576 :     (stored_datum->second)->store_initial_solution();
      76             :   }
      77       16992 :   else if((stored_datum->second)->get_previously_stored() == false) // If we are not at the initial time, we are either creating a new entry or overwriting an existing one
      78             :   {
      79        7980 :     (stored_datum->second)->store_primal_solution(stored_datum);
      80             :   }
      81             :   else // We are overwriting an existing history data
      82             :   {
      83        8540 :     (stored_datum->second)->rewrite_stored_solution();
      84             :   }
      85       17080 : }
      86             : 
      87       23800 : void MemorySolutionHistory::retrieve(bool is_adjoint_solve, Real time)
      88             : {
      89       23800 :   this->find_stored_entry(time, false);
      90             : 
      91             :   // If we are solving the adjoint, the timestep we need to move to the past step
      92             :   // is the one taken at that step to get to the current time.
      93             :   // At the initial time, be ready for the primal time march again.
      94       23800 :   if(is_adjoint_solve)
      95             :   {
      96        8540 :     if( stored_datum != stored_data.begin() )
      97             :     {
      98         212 :       stored_data_iterator stored_datum_past = stored_datum;
      99         212 :       stored_datum_past--;
     100             : 
     101        7420 :       _system.deltat = (stored_datum_past->second)->get_deltat_at();
     102             :     }
     103             :     else
     104             :     {
     105        1120 :       _system.deltat = (stored_datum->second)->get_deltat_at();
     106             :     }
     107             :   }
     108             :   else
     109             :   {
     110       15260 :     if( stored_datum != std::prev(stored_data.end()) )
     111       14280 :       _system.deltat = (stored_datum->second)->get_deltat_at();
     112             :     else
     113             :     {
     114          28 :       stored_data_iterator stored_datum_past = stored_datum;
     115          28 :       stored_datum_past--;
     116             : 
     117         980 :       _system.deltat = (stored_datum_past->second)->get_deltat_at();
     118             :     }
     119             : 
     120             :   }
     121             : 
     122             :   // Get the time at which we are recovering the solution vectors
     123       23800 :   Real recovery_time = stored_datum->first;
     124             : 
     125             :   // Print out what time we are recovering vectors at
     126             :   //    libMesh::out << "Recovering solution vectors at time: " <<
     127             :   //                 recovery_time << std::endl;
     128             : 
     129             :   // Do we not have a solution for this time?  Then
     130             :   // there's nothing to do.
     131       24480 :   if (stored_datum == stored_data.end() ||
     132       23800 :       std::abs(recovery_time - time) > TOLERANCE)
     133             :     {
     134             :       //libMesh::out << "No more solutions to recover ! We are at time t = " <<
     135             :       //                     _system.time << std::endl;
     136          16 :       return;
     137             :     }
     138             : 
     139       23240 :   (stored_datum->second)->retrieve_primal_solution();
     140             : 
     141             :   // We need to call update to put system in a consistent state
     142             :   // with the solution that was read in
     143       23240 :   _system.update();
     144             : }
     145             : 
     146             : }