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         515 : OptimizationFunctionInnerProductHelper::validParams()
      22             : {
      23         515 :   InputParameters params = emptyInputParameters();
      24        1030 :   params.addRequiredParam<FunctionName>("function", "Optimization function.");
      25        1030 :   params.addParam<Real>(
      26             :       "reverse_time_end",
      27        1030 :       0.0,
      28             :       "End time used for reversing the time integration when evaluating function derivative.");
      29         515 :   return params;
      30           0 : }
      31             : 
      32         260 : OptimizationFunctionInnerProductHelper::OptimizationFunctionInnerProductHelper(
      33         260 :     const InputParameters & parameters)
      34         260 :   : _ip_problem(*parameters.getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
      35         260 :     _function(dynamic_cast<const OptimizationFunction *>(&_ip_problem.getFunction(
      36         520 :         parameters.get<FunctionName>("function"), parameters.get<THREAD_ID>("_tid")))),
      37         260 :     _reverse_time_end(parameters.get<Real>("reverse_time_end")),
      38         260 :     _simulation_time(_ip_problem.getMooseApp().getStartTime())
      39             : {
      40         260 :   if (!_function)
      41           0 :     parameters.paramError("function", "Function must be an OptimizationFunction.");
      42         260 : }
      43             : 
      44             : void
      45        7776 : 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        7776 :   _curr_time_ip = nullptr;
      50       55191 :   for (auto & pr : _time_ip)
      51       47415 :     if (MooseUtils::relativeFuzzyEqual(time, pr.first))
      52        7150 :       _curr_time_ip = &pr.second;
      53        7776 :   if (!_curr_time_ip)
      54             :   {
      55         626 :     _time_ip.emplace_back(time, std::vector<Real>());
      56         626 :     _curr_time_ip = &_time_ip.back().second;
      57             :   }
      58        7776 :   _curr_time_ip->clear();
      59             : 
      60        7776 :   _actual_time =
      61        7776 :       MooseUtils::absoluteFuzzyEqual(_reverse_time_end, 0.0) ? time : _reverse_time_end - time + dt;
      62        7776 : }
      63             : 
      64             : void
      65     2618476 : OptimizationFunctionInnerProductHelper::update(const Point & q_point, Real q_inner_product)
      66             : {
      67     2618476 :   if (!_curr_time_ip)
      68           0 :     mooseError("Internal error, 'setCurrentTime' needs to be called before calling 'update'.");
      69             : 
      70     2618476 :   const std::vector<Real> pg = _function->parameterGradient(_actual_time, q_point);
      71     2618476 :   _curr_time_ip->resize(std::max(pg.size(), _curr_time_ip->size()), 0.0);
      72   126163188 :   for (const auto & i : index_range(pg))
      73   123544712 :     (*_curr_time_ip)[i] += q_inner_product * pg[i];
      74     2618476 : }
      75             : 
      76             : void
      77          26 : OptimizationFunctionInnerProductHelper::add(const OptimizationFunctionInnerProductHelper & other)
      78             : {
      79          26 :   _curr_time_ip->resize(std::max(_curr_time_ip->size(), other._curr_time_ip->size()), 0.0);
      80         486 :   for (const auto & i : index_range(*other._curr_time_ip))
      81         460 :     (*_curr_time_ip)[i] += (*other._curr_time_ip)[i];
      82          26 : }
      83             : 
      84             : void
      85        7750 : OptimizationFunctionInnerProductHelper::getVector(std::vector<Real> & result)
      86             : {
      87        7750 :   std::size_t nvar = _curr_time_ip->size();
      88        7750 :   _ip_problem.comm().max(nvar);
      89        7750 :   _curr_time_ip->resize(nvar, 0.0);
      90        7750 :   _ip_problem.comm().sum(*_curr_time_ip);
      91             : 
      92        7750 :   if (!_ip_problem.isTransient())
      93        3093 :     result = (*_curr_time_ip);
      94             :   else
      95             :   {
      96             :     // Make sure everything is the same size
      97       49479 :     for (const auto & it : _time_ip)
      98       44822 :       nvar = std::max(nvar, it.second.size());
      99       49479 :     for (auto & it : _time_ip)
     100       44822 :       it.second.resize(nvar);
     101        4657 :     result.assign(nvar, 0.0);
     102             : 
     103             :     // Integrate in time
     104        4657 :     std::sort(_time_ip.begin(),
     105             :               _time_ip.end(),
     106             :               [](const std::pair<Real, std::vector<Real>> & a,
     107       80096 :                  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       49479 :     for (const auto & ti : _time_ip)
     113      289900 :       for (const auto & i : make_range(nvar))
     114      245078 :         result[i] += ti.second[i];
     115             :   }
     116        7750 : }