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 "OptimizationFunctionInnerProductHelper.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "InputParameters.h"
      14             : #include "MooseError.h"
      15             : #include "FEProblemBase.h"
      16             : 
      17             : // Optimization includes
      18             : #include "OptimizationFunction.h"
      19             : 
      20             : InputParameters
      21         717 : OptimizationFunctionInnerProductHelper::validParams()
      22             : {
      23         717 :   InputParameters params = emptyInputParameters();
      24        1434 :   params.addRequiredParam<FunctionName>("function", "Optimization function.");
      25        1434 :   params.addParam<Real>(
      26             :       "reverse_time_end",
      27        1434 :       0.0,
      28             :       "End time used for reversing the time integration when evaluating function derivative.");
      29         717 :   return params;
      30           0 : }
      31             : 
      32         363 : OptimizationFunctionInnerProductHelper::OptimizationFunctionInnerProductHelper(
      33         363 :     const InputParameters & parameters)
      34         363 :   : _ip_problem(*parameters.getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
      35         363 :     _function(dynamic_cast<const OptimizationFunction *>(&_ip_problem.getFunction(
      36         726 :         parameters.get<FunctionName>("function"), parameters.get<THREAD_ID>("_tid")))),
      37         363 :     _reverse_time_end(parameters.get<Real>("reverse_time_end")),
      38         363 :     _simulation_time(_ip_problem.getMooseApp().getStartTime())
      39             : {
      40         363 :   if (!_function)
      41           0 :     parameters.paramError("function", "Function must be an OptimizationFunction.");
      42         363 : }
      43             : 
      44             : void
      45       10017 : OptimizationFunctionInnerProductHelper::setCurrentTime(Real time, Real dt)
      46             : {
      47             :   // We are saving the integral for each time step.
      48             :   // So create or grab the vector if the time has changed or this is our first time here.
      49       10017 :   _curr_time_ip = nullptr;
      50       79653 :   for (auto & pr : _time_ip)
      51       69636 :     if (MooseUtils::relativeFuzzyEqual(time, pr.first))
      52        9081 :       _curr_time_ip = &pr.second;
      53       10017 :   if (!_curr_time_ip)
      54             :   {
      55         936 :     _time_ip.emplace_back(time, std::vector<Real>());
      56         936 :     _curr_time_ip = &_time_ip.back().second;
      57             :   }
      58       10017 :   _curr_time_ip->clear();
      59             : 
      60       10017 :   _actual_time =
      61       10017 :       MooseUtils::absoluteFuzzyEqual(_reverse_time_end, 0.0) ? time : _reverse_time_end - time + dt;
      62       10017 : }
      63             : 
      64             : void
      65     2948378 : OptimizationFunctionInnerProductHelper::update(const Point & q_point, Real q_inner_product)
      66             : {
      67     2948378 :   if (!_curr_time_ip)
      68           0 :     mooseError("Internal error, 'setCurrentTime' needs to be called before calling 'update'.");
      69             : 
      70     2948378 :   const std::vector<Real> pg = _function->parameterGradient(_actual_time, q_point);
      71     2948378 :   _curr_time_ip->resize(std::max(pg.size(), _curr_time_ip->size()), 0.0);
      72    20902222 :   for (const auto & i : index_range(pg))
      73    17953844 :     (*_curr_time_ip)[i] += q_inner_product * pg[i];
      74     2948378 : }
      75             : 
      76             : void
      77          66 : OptimizationFunctionInnerProductHelper::add(const OptimizationFunctionInnerProductHelper & other)
      78             : {
      79          66 :   _curr_time_ip->resize(std::max(_curr_time_ip->size(), other._curr_time_ip->size()), 0.0);
      80        1434 :   for (const auto & i : index_range(*other._curr_time_ip))
      81        1368 :     (*_curr_time_ip)[i] += (*other._curr_time_ip)[i];
      82          66 : }
      83             : 
      84             : void
      85        9951 : OptimizationFunctionInnerProductHelper::getVector(std::vector<Real> & result)
      86             : {
      87        9951 :   std::size_t nvar = _curr_time_ip->size();
      88        9951 :   _ip_problem.comm().max(nvar);
      89        9951 :   _curr_time_ip->resize(nvar, 0.0);
      90        9951 :   _ip_problem.comm().sum(*_curr_time_ip);
      91             : 
      92        9951 :   if (!_ip_problem.isTransient())
      93        2962 :     result = (*_curr_time_ip);
      94             :   else
      95             :   {
      96             :     // Make sure everything is the same size
      97       74233 :     for (const auto & it : _time_ip)
      98       67244 :       nvar = std::max(nvar, it.second.size());
      99       74233 :     for (auto & it : _time_ip)
     100       67244 :       it.second.resize(nvar);
     101        6989 :     result.assign(nvar, 0.0);
     102             : 
     103             :     // Integrate in time
     104        6989 :     std::sort(_time_ip.begin(),
     105             :               _time_ip.end(),
     106             :               [](const std::pair<Real, std::vector<Real>> & a,
     107      120141 :                  const std::pair<Real, std::vector<Real>> & b) { return a.first < b.first; });
     108             :     // We are integrating over the entire history here. Technically this can be done in an
     109             :     // accumulative manner by storing the integration over previous time steps. However, this
     110             :     // is a relatively cheap operation and considering all cases where time steps may be
     111             :     // repeated would be difficult.
     112       74233 :     for (const auto & ti : _time_ip)
     113      434575 :       for (const auto & i : make_range(nvar))
     114      367331 :         result[i] += ti.second[i];
     115             :   }
     116        9951 : }