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 "ParameterMeshFunction.h"
      11             : 
      12             : #include "AddVariableAction.h"
      13             : 
      14             : registerMooseObject("OptimizationApp", ParameterMeshFunction);
      15             : 
      16             : InputParameters
      17         380 : ParameterMeshFunction::validParams()
      18             : {
      19         380 :   InputParameters params = OptimizationFunction::validParams();
      20         380 :   params.addClassDescription("Optimization function with parameters represented by a mesh and "
      21             :                              "finite-element shape functions.");
      22             : 
      23         760 :   params.addRequiredParam<FileName>("exodus_mesh", "File containing parameter mesh.");
      24         760 :   params.addParam<MooseEnum>("family",
      25         760 :                              AddVariableAction::getNonlinearVariableFamilies(),
      26             :                              "Family of FE shape functions for parameter.");
      27         760 :   params.addParam<MooseEnum>("order",
      28         760 :                              AddVariableAction::getNonlinearVariableOrders(),
      29             :                              "Order of FE shape functions for parameter.");
      30             : 
      31         760 :   params.addRequiredParam<ReporterName>(
      32             :       "parameter_name", "Reporter or VectorPostprocessor vector containing parameter values.");
      33         760 :   params.addParam<ReporterName>("time_name",
      34             :                                 "Name of vector-postprocessor or reporter vector containing time, "
      35             :                                 "default assumes time independence.");
      36         760 :   params.addParam<bool>(
      37         760 :       "project_points", false, "Whether to find the closest point on parameter mesh.");
      38         760 :   params.addParam<unsigned int>(
      39             :       "kdtree_candidates",
      40         760 :       5,
      41             :       "Number of nearest node candidates to consider when projecting points to parameter mesh.");
      42             : 
      43         380 :   return params;
      44           0 : }
      45             : 
      46         196 : ParameterMeshFunction::ParameterMeshFunction(const InputParameters & parameters)
      47             :   : OptimizationFunction(parameters),
      48             :     ReporterInterface(this),
      49         196 :     _parameter_mesh(AddVariableAction::feType(parameters),
      50         194 :                     getParam<FileName>("exodus_mesh"),
      51         586 :                     getParam<bool>("project_points"),
      52         390 :                     getParam<unsigned int>("kdtree_candidates")),
      53         194 :     _values(getReporterValue<std::vector<Real>>("parameter_name")),
      54         398 :     _coordt(isParamValid("time_name") ? getReporterValue<std::vector<Real>>("time_name")
      55         196 :                                       : _empty_vec)
      56             : {
      57         194 : }
      58             : 
      59             : Real
      60     7606510 : ParameterMeshFunction::value(Real t, const Point & p) const
      61             : {
      62     7606510 :   checkSize();
      63             : 
      64     7606509 :   const auto ti = interpolateTime(t);
      65     7606509 :   const dof_id_type offset0 = ti[0].first * _parameter_mesh.size();
      66     7606509 :   const dof_id_type offset1 = ti[1].first * _parameter_mesh.size();
      67             : 
      68             :   std::vector<dof_id_type> dof_indices;
      69             :   std::vector<Real> weights;
      70             : 
      71     7606509 :   _parameter_mesh.getIndexAndWeight(p, dof_indices, weights);
      72             : 
      73             :   Real val = 0;
      74    38394739 :   for (const auto & i : index_range(dof_indices))
      75    30788230 :     val += (_values[dof_indices[i] + offset0] * ti[0].second +
      76    30788230 :             _values[dof_indices[i] + offset1] * ti[1].second) *
      77             :            weights[i];
      78     7606509 :   return val;
      79     7606509 : }
      80             : 
      81             : RealGradient
      82       10800 : ParameterMeshFunction::gradient(Real t, const Point & p) const
      83             : {
      84       10800 :   checkSize();
      85             : 
      86       10800 :   const auto ti = interpolateTime(t);
      87       10800 :   const dof_id_type offset0 = ti[0].first * _parameter_mesh.size();
      88       10800 :   const dof_id_type offset1 = ti[1].first * _parameter_mesh.size();
      89             : 
      90             :   std::vector<dof_id_type> dof_indices;
      91             :   std::vector<RealGradient> weights;
      92       10800 :   _parameter_mesh.getIndexAndWeight(p, dof_indices, weights);
      93             : 
      94             :   RealGradient val(0, 0, 0);
      95       56400 :   for (const auto & i : index_range(dof_indices))
      96       45600 :     val += (_values[dof_indices[i] + offset0] * ti[0].second +
      97       45600 :             _values[dof_indices[i] + offset1] * ti[1].second) *
      98             :            weights[i];
      99       10800 :   return val;
     100       10800 : }
     101             : 
     102             : Real
     103           0 : ParameterMeshFunction::timeDerivative(Real t, const Point & p) const
     104             : {
     105           0 :   checkSize();
     106             : 
     107           0 :   const auto ti = interpolateTime(t, true);
     108           0 :   if (ti[0].first == ti[1].first)
     109             :     return 0.0;
     110           0 :   const dof_id_type offset0 = ti[0].first * _parameter_mesh.size();
     111           0 :   const dof_id_type offset1 = ti[1].first * _parameter_mesh.size();
     112             : 
     113             :   std::vector<dof_id_type> dof_indices;
     114             :   std::vector<Real> weights;
     115           0 :   _parameter_mesh.getIndexAndWeight(p, dof_indices, weights);
     116             : 
     117             :   Real val = 0;
     118           0 :   for (const auto & i : index_range(dof_indices))
     119           0 :     val += (_values[dof_indices[i] + offset0] * ti[0].second +
     120           0 :             _values[dof_indices[i] + offset1] * ti[1].second) *
     121             :            weights[i];
     122             :   return val;
     123           0 : }
     124             : 
     125             : std::vector<Real>
     126     1046134 : ParameterMeshFunction::parameterGradient(Real t, const Point & p) const
     127             : {
     128     1046134 :   const auto ti = interpolateTime(t);
     129     1046134 :   const dof_id_type offset0 = ti[0].first * _parameter_mesh.size();
     130     1046134 :   const dof_id_type offset1 = ti[1].first * _parameter_mesh.size();
     131             : 
     132             :   std::vector<dof_id_type> dof_indices;
     133             :   std::vector<Real> weights;
     134     1046134 :   _parameter_mesh.getIndexAndWeight(p, dof_indices, weights);
     135             : 
     136             :   dof_id_type sz = _parameter_mesh.size();
     137     1046134 :   if (!_coordt.empty())
     138        4356 :     sz *= _coordt.size();
     139     1046134 :   std::vector<Real> pg(sz, 0.0);
     140     5270922 :   for (const auto & i : index_range(dof_indices))
     141             :   {
     142     4224788 :     pg[dof_indices[i] + offset0] += weights[i] * ti[0].second;
     143     4224788 :     pg[dof_indices[i] + offset1] += weights[i] * ti[1].second;
     144             :   }
     145     1046134 :   return pg;
     146     1046134 : }
     147             : 
     148             : std::array<std::pair<std::size_t, Real>, 2>
     149     8663443 : ParameterMeshFunction::interpolateTime(Real t, bool derivative) const
     150             : {
     151     8663443 :   std::array<std::pair<std::size_t, Real>, 2> ti;
     152     8663443 :   if (_coordt.size() <= 1 || MooseUtils::absoluteFuzzyLessEqual(t, _coordt[0]))
     153             :   {
     154             :     ti[0] = {0, 0.5};
     155             :     ti[1] = {0, 0.5};
     156             :   }
     157       12630 :   else if (MooseUtils::absoluteFuzzyGreaterEqual(t, _coordt.back()))
     158             :   {
     159        2526 :     ti[0] = {_coordt.size() - 1, 0.5};
     160             :     ti[1] = {_coordt.size() - 1, 0.5};
     161             :   }
     162             :   else
     163       15156 :     for (std::size_t i = 1; i < _coordt.size(); ++i)
     164       15156 :       if (MooseUtils::absoluteFuzzyGreaterEqual(_coordt[i], t))
     165             :       {
     166       10104 :         const Real dt = _coordt[i] - _coordt[i - 1];
     167       10104 :         if (derivative)
     168             :         {
     169           0 :           ti[0] = {i - 1, -1.0 / dt};
     170             :           ti[1] = {i, 1.0 / dt};
     171             :         }
     172             :         else
     173             :         {
     174       10104 :           ti[0] = {i - 1, (_coordt[i] - t) / dt};
     175       10104 :           ti[1] = {i, (t - _coordt[i - 1]) / dt};
     176             :         }
     177             :         break;
     178             :       }
     179             : 
     180     8663443 :   return ti;
     181             : }
     182             : 
     183             : void
     184     7617310 : ParameterMeshFunction::checkSize() const
     185             : {
     186             :   dof_id_type sz = _parameter_mesh.size();
     187     7617310 :   if (!_coordt.empty())
     188       10800 :     sz *= _coordt.size();
     189     7617310 :   if (sz != _values.size())
     190           1 :     paramError("parameter_name",
     191             :                "Size of parameter vector (",
     192             :                _values.size(),
     193             :                ") does not match number of degrees of freedom in mesh (",
     194             :                sz,
     195             :                ").");
     196     7617309 : }