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             : 
      19             : #include "libmesh/diff_context.h"
      20             : #include "libmesh/diff_system.h"
      21             : #include "libmesh/diff_system.h"
      22             : #include "libmesh/unsteady_solver.h"
      23             : 
      24             : // C++ includes
      25             : #include <memory>
      26             : 
      27             : namespace libMesh
      28             : {
      29             : 
      30             : 
      31             : 
      32     2011429 : DiffContext::DiffContext (const System & sys,
      33     2011429 :                           bool allocate_local_matrices) :
      34     2011429 :   time(sys.time),
      35     1928550 :   system_time(sys.time),
      36     1845671 :   elem_solution_derivative(1.),
      37     1845671 :   elem_solution_rate_derivative(1.),
      38     1845671 :   elem_solution_accel_derivative(1.),
      39     1845671 :   fixed_solution_derivative(0.),
      40     1845671 :   _have_local_matrices(allocate_local_matrices),
      41     1928550 :   _dof_indices_var(sys.n_vars()),
      42     1845671 :   _deltat(nullptr),
      43     1845671 :   _system(sys),
      44     2177187 :   _is_adjoint(false)
      45             : {
      46             :   // Finally initialize solution/residual/jacobian data structures
      47       82879 :   unsigned int nv = sys.n_vars();
      48             : 
      49     2011429 :   _elem_subsolutions.reserve(nv);
      50     2011429 :   _elem_subresiduals.reserve(nv);
      51     2011429 :   _elem_subsolution_rates.reserve(nv);
      52     2011429 :   _elem_subsolution_accels.reserve(nv);
      53             : 
      54     2011429 :   if (sys.use_fixed_solution)
      55           0 :     _elem_fixed_subsolutions.reserve(nv);
      56             : 
      57     2011429 :   if (allocate_local_matrices)
      58      350510 :     _elem_subjacobians.resize(nv);
      59             : 
      60             :   // If the user resizes sys.qoi, it will invalidate us
      61     2011429 :   std::size_t n_qoi = sys.n_qois();
      62     2011429 :   _elem_qoi.resize(n_qoi);
      63     2011429 :   _elem_qoi_derivative.resize(n_qoi);
      64     2011429 :   _elem_qoi_subderivatives.resize(n_qoi);
      65     2584116 :   for (std::size_t q=0; q != n_qoi; ++q)
      66      605971 :     _elem_qoi_subderivatives[q].reserve(nv);
      67             : 
      68     4407969 :   for (unsigned int i=0; i != nv; ++i)
      69             :     {
      70     2396540 :       _elem_subsolutions.emplace_back(_elem_solution);
      71     2396540 :       _elem_subresiduals.emplace_back(_elem_residual);
      72     3052315 :       for (std::size_t q=0; q != n_qoi; ++q)
      73      733911 :         _elem_qoi_subderivatives[q].emplace_back(_elem_qoi_derivative[q]);
      74     2396540 :       if (allocate_local_matrices)
      75      556948 :         _elem_subjacobians[i].reserve(nv);
      76             : 
      77             :       // Only make space for these if we're using DiffSystem
      78             :       // This is assuming *only* DiffSystem is using elem_solution_rate/accel
      79     2396540 :       const DifferentiableSystem * diff_system = dynamic_cast<const DifferentiableSystem *>(&sys);
      80     2396540 :       if (diff_system)
      81             :         {
      82             :           // Now, we only need these if the solver is unsteady
      83      710516 :           if (!diff_system->get_time_solver().is_steady())
      84             :             {
      85      387720 :               _elem_subsolution_rates.emplace_back(_elem_solution_rate);
      86             : 
      87             :               // We only need accel space if the TimeSolver is second order
      88       18430 :               const UnsteadySolver & time_solver = cast_ref<const UnsteadySolver &>(diff_system->get_time_solver());
      89             : 
      90      387720 :               if (time_solver.time_order() >= 2 || !diff_system->get_second_order_vars().empty())
      91       43899 :                 _elem_subsolution_accels.emplace_back(_elem_solution_accel);
      92             :             }
      93             :         }
      94             : 
      95     2396540 :       if (sys.use_fixed_solution)
      96           0 :         _elem_fixed_subsolutions.emplace_back(_elem_fixed_solution);
      97             : 
      98     2396540 :       if (allocate_local_matrices)
      99     1689224 :         for (unsigned int j=0; j != nv; ++j)
     100     1208576 :           _elem_subjacobians[i].emplace_back(_elem_jacobian);
     101             :     }
     102     2011429 : }
     103             : 
     104             : 
     105             : 
     106    12919697 : DiffContext::~DiffContext () = default;
     107             : 
     108             : 
     109      492769 : void DiffContext::set_deltat_pointer(Real * dt)
     110             : {
     111             :   // We may actually want to be able to set this pointer to nullptr, so
     112             :   // don't report an error for that.
     113      492769 :   _deltat = dt;
     114      492769 : }
     115             : 
     116             : 
     117    43285856 : Real DiffContext::get_deltat_value()
     118             : {
     119     3887976 :   libmesh_assert(_deltat);
     120             : 
     121    43285856 :   return *_deltat;
     122             : }
     123             : 
     124             : 
     125      186048 : void DiffContext::add_localized_vector (NumericVector<Number> & localized_vector, const System & sys)
     126             : {
     127             :   // Make an empty pair keyed with a reference to this _localized_vector
     128      362304 :   _localized_vectors[&localized_vector] = std::make_pair(DenseVector<Number>(), std::vector<DenseSubVector<Number>>());
     129             : 
     130        9792 :   unsigned int nv = sys.n_vars();
     131             : 
     132      186048 :   _localized_vectors[&localized_vector].second.reserve(nv);
     133             : 
     134             :   // Fill the DenseSubVector with nv copies of DenseVector
     135      372096 :   for (unsigned int i=0; i != nv; ++i)
     136      186048 :     _localized_vectors[&localized_vector].second.emplace_back(_localized_vectors[&localized_vector].first);
     137      186048 : }
     138             : 
     139             : 
     140           0 : DenseVector<Number> & DiffContext::get_localized_vector (const NumericVector<Number> & localized_vector)
     141             : {
     142           0 :   return _localized_vectors[&localized_vector].first;
     143             : }
     144             : 
     145             : 
     146           0 : const DenseVector<Number> & DiffContext::get_localized_vector (const NumericVector<Number> & localized_vector) const
     147             : {
     148           0 :   auto localized_vectors_it = _localized_vectors.find(&localized_vector);
     149           0 :   libmesh_assert(localized_vectors_it != _localized_vectors.end());
     150           0 :   return localized_vectors_it->second.first;
     151             : }
     152             : 
     153             : 
     154           0 : DenseSubVector<Number> & DiffContext::get_localized_subvector (const NumericVector<Number> & localized_vector, unsigned int var)
     155             : {
     156           0 :   return _localized_vectors[&localized_vector].second[var];
     157             : }
     158             : 
     159             : 
     160     2173248 : const DenseSubVector<Number> & DiffContext::get_localized_subvector (const NumericVector<Number> & localized_vector, unsigned int var) const
     161             : {
     162     2173248 :   auto localized_vectors_it = _localized_vectors.find(&localized_vector);
     163      197568 :   libmesh_assert(localized_vectors_it != _localized_vectors.end());
     164     2568384 :   return localized_vectors_it->second.second[var];
     165             : }
     166             : 
     167             : } // namespace libMesh