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         680 : ParameterMeshFunction::validParams()
      18             : {
      19         680 :   InputParameters params = OptimizationFunction::validParams();
      20         680 :   params.addClassDescription("Optimization function with parameters represented by a mesh and "
      21             :                              "finite-element shape functions.");
      22             : 
      23        1360 :   params.addRequiredParam<FileName>("exodus_mesh", "File containing parameter mesh.");
      24        1360 :   params.addParam<MooseEnum>("family",
      25        1360 :                              AddVariableAction::getNonlinearVariableFamilies(),
      26             :                              "Family of FE shape functions for parameter.");
      27        1360 :   params.addParam<MooseEnum>("order",
      28        1360 :                              AddVariableAction::getNonlinearVariableOrders(),
      29             :                              "Order of FE shape functions for parameter.");
      30             : 
      31        1360 :   params.addRequiredParam<ReporterName>(
      32             :       "parameter_name", "Reporter or VectorPostprocessor vector containing parameter values.");
      33        1360 :   params.addParam<ReporterName>("time_name",
      34             :                                 "Name of vector-postprocessor or reporter vector containing time, "
      35             :                                 "default assumes time independence.");
      36        1360 :   params.addParam<bool>(
      37        1360 :       "project_points", false, "Whether to find the closest point on parameter mesh.");
      38        1360 :   params.addParam<unsigned int>(
      39             :       "kdtree_candidates",
      40        1360 :       5,
      41             :       "Number of nearest node candidates to consider when projecting points to parameter mesh.");
      42             : 
      43         680 :   return params;
      44           0 : }
      45             : 
      46         356 : ParameterMeshFunction::ParameterMeshFunction(const InputParameters & parameters)
      47             :   : OptimizationFunction(parameters),
      48             :     ReporterInterface(this),
      49         356 :     _parameter_mesh(AddVariableAction::feType(parameters),
      50         354 :                     getParam<FileName>("exodus_mesh"),
      51        1066 :                     getParam<bool>("project_points"),
      52         710 :                     getParam<unsigned int>("kdtree_candidates")),
      53         354 :     _values(getReporterValue<std::vector<Real>>("parameter_name")),
      54         730 :     _coordt(isParamValid("time_name") ? getReporterValue<std::vector<Real>>("time_name")
      55         356 :                                       : _empty_vec)
      56             : {
      57         354 : }
      58             : 
      59             : Real
      60    11405013 : ParameterMeshFunction::value(Real t, const Point & p) const
      61             : {
      62    11405013 :   checkSize();
      63             : 
      64    11405011 :   const auto ti = interpolateTime(t);
      65    11405011 :   const dof_id_type offset0 = ti[0].first * _parameter_mesh.size();
      66    11405011 :   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    11405011 :   _parameter_mesh.getIndexAndWeight(p, dof_indices, weights);
      72             : 
      73             :   Real val = 0;
      74    57560360 :   for (const auto & i : index_range(dof_indices))
      75    46155349 :     val += (_values[dof_indices[i] + offset0] * ti[0].second +
      76    46155349 :             _values[dof_indices[i] + offset1] * ti[1].second) *
      77             :            weights[i];
      78    11405011 :   return val;
      79    11405011 : }
      80             : 
      81             : RealGradient
      82       16200 : ParameterMeshFunction::gradient(Real t, const Point & p) const
      83             : {
      84       16200 :   checkSize();
      85             : 
      86       16200 :   const auto ti = interpolateTime(t);
      87       16200 :   const dof_id_type offset0 = ti[0].first * _parameter_mesh.size();
      88       16200 :   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       16200 :   _parameter_mesh.getIndexAndWeight(p, dof_indices, weights);
      93             : 
      94             :   RealGradient val(0, 0, 0);
      95       84600 :   for (const auto & i : index_range(dof_indices))
      96       68400 :     val += (_values[dof_indices[i] + offset0] * ti[0].second +
      97       68400 :             _values[dof_indices[i] + offset1] * ti[1].second) *
      98             :            weights[i];
      99       16200 :   return val;
     100       16200 : }
     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     1569401 : ParameterMeshFunction::parameterGradient(Real t, const Point & p) const
     127             : {
     128     1569401 :   const auto ti = interpolateTime(t);
     129     1569401 :   const dof_id_type offset0 = ti[0].first * _parameter_mesh.size();
     130     1569401 :   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     1569401 :   _parameter_mesh.getIndexAndWeight(p, dof_indices, weights);
     135             : 
     136             :   dof_id_type sz = _parameter_mesh.size();
     137     1569401 :   if (!_coordt.empty())
     138        6534 :     sz *= _coordt.size();
     139     1569401 :   std::vector<Real> pg(sz, 0.0);
     140     7908383 :   for (const auto & i : index_range(dof_indices))
     141             :   {
     142     6338982 :     pg[dof_indices[i] + offset0] += weights[i] * ti[0].second;
     143     6338982 :     pg[dof_indices[i] + offset1] += weights[i] * ti[1].second;
     144             :   }
     145     1569401 :   return pg;
     146     1569401 : }
     147             : 
     148             : std::array<std::pair<std::size_t, Real>, 2>
     149    12990612 : ParameterMeshFunction::interpolateTime(Real t, bool derivative) const
     150             : {
     151    12990612 :   std::array<std::pair<std::size_t, Real>, 2> ti;
     152    12990612 :   if (_coordt.size() <= 1 || MooseUtils::absoluteFuzzyLessEqual(t, _coordt[0]))
     153             :   {
     154             :     ti[0] = {0, 0.5};
     155             :     ti[1] = {0, 0.5};
     156             :   }
     157       18945 :   else if (MooseUtils::absoluteFuzzyGreaterEqual(t, _coordt.back()))
     158             :   {
     159        3789 :     ti[0] = {_coordt.size() - 1, 0.5};
     160             :     ti[1] = {_coordt.size() - 1, 0.5};
     161             :   }
     162             :   else
     163       22734 :     for (std::size_t i = 1; i < _coordt.size(); ++i)
     164       22734 :       if (MooseUtils::absoluteFuzzyGreaterEqual(_coordt[i], t))
     165             :       {
     166       15156 :         const Real dt = _coordt[i] - _coordt[i - 1];
     167       15156 :         if (derivative)
     168             :         {
     169           0 :           ti[0] = {i - 1, -1.0 / dt};
     170             :           ti[1] = {i, 1.0 / dt};
     171             :         }
     172             :         else
     173             :         {
     174       15156 :           ti[0] = {i - 1, (_coordt[i] - t) / dt};
     175       15156 :           ti[1] = {i, (t - _coordt[i - 1]) / dt};
     176             :         }
     177             :         break;
     178             :       }
     179             : 
     180    12990612 :   return ti;
     181             : }
     182             : 
     183             : void
     184    11421213 : ParameterMeshFunction::checkSize() const
     185             : {
     186             :   dof_id_type sz = _parameter_mesh.size();
     187    11421213 :   if (!_coordt.empty())
     188       16200 :     sz *= _coordt.size();
     189    11421213 :   if (sz != _values.size())
     190           2 :     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    11421211 : }