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_solver.h"
      20             : #include "libmesh/diff_system.h"
      21             : #include "libmesh/steady_solver.h"
      22             : 
      23             : #include "libmesh/adjoint_refinement_estimator.h"
      24             : #include "libmesh/error_vector.h"
      25             : 
      26             : namespace libMesh
      27             : {
      28             : 
      29             : 
      30             : 
      31        5226 : SteadySolver::~SteadySolver () = default;
      32             : 
      33             : 
      34             : 
      35     5296989 : bool SteadySolver::element_residual(bool request_jacobian,
      36             :                                     DiffContext & context)
      37             : {
      38     5789111 :   return this->_general_residual(request_jacobian,
      39             :                                  context,
      40             :                                  &DifferentiablePhysics::element_time_derivative,
      41     5296989 :                                  &DifferentiablePhysics::element_constraint);
      42             : }
      43             : 
      44             : 
      45             : 
      46     1082026 : bool SteadySolver::side_residual(bool request_jacobian,
      47             :                                  DiffContext & context)
      48             : {
      49     1200592 :   return this->_general_residual(request_jacobian,
      50             :                                  context,
      51             :                                  &DifferentiablePhysics::side_time_derivative,
      52     1082026 :                                  &DifferentiablePhysics::side_constraint);
      53             : }
      54             : 
      55             : 
      56             : 
      57           0 : bool SteadySolver::nonlocal_residual(bool request_jacobian,
      58             :                                      DiffContext & context)
      59             : {
      60           0 :   return this->_general_residual(request_jacobian,
      61             :                                  context,
      62             :                                  &DifferentiablePhysics::nonlocal_time_derivative,
      63           0 :                                  &DifferentiablePhysics::nonlocal_constraint);
      64             : }
      65             : 
      66             : 
      67             : 
      68     6379015 : bool SteadySolver::_general_residual(bool request_jacobian,
      69             :                                      DiffContext & context,
      70             :                                      ResFuncType time_deriv,
      71             :                                      ResFuncType constraint)
      72             : {
      73             :   // If a fixed solution is requested, it will just be the current
      74             :   // solution
      75     6379015 :   if (_system.use_fixed_solution)
      76             :     {
      77           0 :       context.get_elem_fixed_solution() = context.get_elem_solution();
      78           0 :       context.fixed_solution_derivative = 1.0;
      79             :     }
      80             : 
      81             :   bool jacobian_computed =
      82    12147342 :     (_system.get_physics()->*time_deriv)(request_jacobian, context);
      83             : 
      84             :   // The user shouldn't compute a jacobian unless requested
      85      610688 :   libmesh_assert (request_jacobian || !jacobian_computed);
      86             : 
      87             :   bool jacobian_computed2 =
      88    12147342 :     (_system.get_physics()->*constraint)(jacobian_computed, context);
      89             : 
      90             :   // The user shouldn't compute a jacobian unless requested
      91      610688 :   libmesh_assert (jacobian_computed || !jacobian_computed2);
      92             : 
      93     6379015 :   return jacobian_computed2;
      94             : }
      95             : 
      96           0 : void SteadySolver::integrate_qoi_timestep()
      97             : {
      98           0 :   this->_system.assemble_qoi();
      99             : 
     100           0 :   return;
     101             : }
     102             : 
     103           0 : void SteadySolver::integrate_adjoint_sensitivity(const QoISet & qois, const ParameterVector & parameter_vector, SensitivityData & sensitivities)
     104             : {
     105           0 :   this->_system.ImplicitSystem::adjoint_qoi_parameter_sensitivity(qois, parameter_vector, sensitivities);
     106             : 
     107           0 :   return;
     108             : }
     109             : 
     110             : #ifdef LIBMESH_ENABLE_AMR
     111           0 : void SteadySolver::integrate_adjoint_refinement_error_estimate
     112             :   (AdjointRefinementEstimator & adjoint_refinement_error_estimator,
     113             :    ErrorVector & QoI_elementwise_error)
     114             : {
     115             :   // Base class assumes a direct steady state error estimate
     116           0 :   adjoint_refinement_error_estimator.estimate_error(_system, QoI_elementwise_error);
     117             : 
     118             :   // Also get the spatially integrated errors for all the QoIs in the QoI set
     119           0 :   for (auto j : make_range(_system.n_qois()))
     120             :   {
     121             :     // Skip this QoI if not in the QoI Set
     122           0 :     if (adjoint_refinement_error_estimator.qoi_set().has_index(j))
     123             :     {
     124           0 :       _system.set_qoi_error_estimate(j, adjoint_refinement_error_estimator.get_global_QoI_error_estimate(j));
     125             :     }
     126             :   }
     127             : 
     128           0 :   return;
     129             : }
     130             : #endif // LIBMESH_ENABLE_AMR
     131             : 
     132             : } // namespace libMesh