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 "FindValueOnLine.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "MooseMesh.h"
      14             : #include "MooseUtils.h"
      15             : #include "MooseVariable.h"
      16             : 
      17             : registerMooseObject("MooseApp", FindValueOnLine);
      18             : 
      19             : InputParameters
      20       14427 : FindValueOnLine::validParams()
      21             : {
      22       14427 :   InputParameters params = GeneralPostprocessor::validParams();
      23       14427 :   params.addClassDescription("Find a specific target value along a sampling line. The variable "
      24             :                              "values along the line should change monotonically. The target value "
      25             :                              "is searched using a bisection algorithm.");
      26       14427 :   params.addParam<Point>("start_point", "Start point of the sampling line.");
      27       14427 :   params.addParam<Point>("end_point", "End point of the sampling line.");
      28       14427 :   params.addParam<Real>("target", "Target value to locate.");
      29       43281 :   params.addParam<bool>(
      30             :       "error_if_not_found",
      31       28854 :       true,
      32             :       "If true, stop with error if target value is not found on the line. If false, "
      33             :       "return default_value.");
      34       43281 :   params.addParam<Real>("default_value",
      35       28854 :                         -1,
      36             :                         "Value to return if target value is not found on line and "
      37             :                         "error_if_not_found is false.");
      38       14427 :   params.addParam<unsigned int>("depth", 36, "Maximum number of bisections to perform.");
      39       43281 :   params.addParam<Real>(
      40             :       "tol",
      41       28854 :       1e-10,
      42             :       "Stop search if a value is found that is equal to the target with this tolerance applied.");
      43       14427 :   params.addCoupledVar("v", "Variable to inspect");
      44       14427 :   return params;
      45           0 : }
      46             : 
      47          81 : FindValueOnLine::FindValueOnLine(const InputParameters & parameters)
      48             :   : GeneralPostprocessor(parameters),
      49             :     Coupleable(this, false),
      50          81 :     _start_point(getParam<Point>("start_point")),
      51          81 :     _end_point(getParam<Point>("end_point")),
      52          81 :     _length((_end_point - _start_point).norm()),
      53          81 :     _target(getParam<Real>("target")),
      54          81 :     _error_if_not_found(getParam<bool>("error_if_not_found")),
      55          81 :     _default_value(getParam<Real>("default_value")),
      56          81 :     _depth(getParam<unsigned int>("depth")),
      57          81 :     _tol(getParam<Real>("tol")),
      58          81 :     _coupled_var(*getVar("v", 0)),
      59          81 :     _position(0.0),
      60          81 :     _mesh(_subproblem.mesh()),
      61         162 :     _point_vec(1)
      62             : {
      63          81 : }
      64             : 
      65             : void
      66         381 : FindValueOnLine::initialize()
      67             : {
      68             :   // We do this here just in case it's been destroyed and recreated becaue of mesh adaptivity.
      69         381 :   _pl = _mesh.getPointLocator();
      70         381 :   _pl->enable_out_of_mesh_mode();
      71         381 : }
      72             : 
      73             : void
      74         381 : FindValueOnLine::execute()
      75             : {
      76             :   Real s;
      77         381 :   Real s_left = 0.0;
      78         381 :   Real left = getValueAtPoint(_start_point);
      79         381 :   Real s_right = 1.0;
      80         381 :   Real right = getValueAtPoint(_end_point);
      81             : 
      82             :   // Helper for erroring; won't error if error_if_not_found=false
      83         364 :   const auto error = [this](auto &&... message)
      84             :   {
      85         182 :     if (_error_if_not_found)
      86          12 :       mooseError(std::scientific, std::setprecision(4), message...);
      87         170 :     _position = _default_value;
      88         547 :   };
      89             : 
      90             :   /**
      91             :    * Here we determine the direction of the solution. i.e. the left might be the high value
      92             :    * while the right might be the low value.
      93             :    */
      94         377 :   bool left_to_right = left < right;
      95             :   // Initial bounds check
      96         377 :   if ((left_to_right && _target < left) || (!left_to_right && _target < right))
      97             :   {
      98         150 :     error("Target value ",
      99         150 :           _target,
     100             :           " is less than the minimum sampled value ",
     101             :           std::min(left, right));
     102         170 :     return;
     103             :   }
     104         227 :   if ((left_to_right && _target > right) || (!left_to_right && _target > left))
     105             :   {
     106          28 :     error("Target value ",
     107          28 :           _target,
     108             :           " is greater than the maximum sampled value ",
     109             :           std::max(left, right));
     110          24 :     return;
     111             :   }
     112             : 
     113         199 :   bool found_it = false;
     114         199 :   Real value = 0;
     115        6612 :   for (unsigned int i = 0; i < _depth; ++i)
     116             :   {
     117             :     // find midpoint
     118        6608 :     s = (s_left + s_right) / 2.0;
     119        6608 :     Point p = s * (_end_point - _start_point) + _start_point;
     120             : 
     121             :     // sample value
     122        6608 :     value = getValueAtPoint(p);
     123             : 
     124             :     // have we hit the target value yet?
     125        6608 :     if (MooseUtils::absoluteFuzzyEqual(value, _target, _tol))
     126             :     {
     127         195 :       found_it = true;
     128         195 :       break;
     129             :     }
     130             : 
     131             :     // bisect
     132        6413 :     if ((left_to_right && _target < value) || (!left_to_right && _target > value))
     133             :       // to the left
     134        2981 :       s_right = s;
     135             :     else
     136             :       // to the right
     137        3432 :       s_left = s;
     138             :   }
     139             : 
     140             :   // Return error if target value (within tol) was not found within depth bisections
     141         199 :   if (!found_it)
     142             :   {
     143           4 :     error("Target value ",
     144           4 :           _target,
     145             :           " not found on line within tolerance.\nLast sample: ",
     146             :           value,
     147             :           ", difference from target: ",
     148           4 :           std::abs(_target - value));
     149           0 :     return;
     150             :   }
     151             : 
     152         195 :   _position = s * _length;
     153             : }
     154             : 
     155             : Real
     156        7370 : FindValueOnLine::getValueAtPoint(const Point & p)
     157             : {
     158        7370 :   const Elem * elem = (*_pl)(p);
     159             : 
     160             :   processor_id_type elem_proc_id =
     161        7370 :       elem ? elem->processor_id() : libMesh::DofObject::invalid_processor_id;
     162        7370 :   _communicator.min(elem_proc_id);
     163             : 
     164        7370 :   if (elem_proc_id == libMesh::DofObject::invalid_processor_id)
     165             :   {
     166             :     // there is no element
     167           4 :     mooseError("No element found at the current search point. Please make sure the sampling line "
     168             :                "stays inside the mesh completely.");
     169             :   }
     170             : 
     171        7366 :   Real value = 0;
     172             : 
     173        7366 :   if (elem)
     174             :   {
     175        6815 :     if (elem->processor_id() == processor_id())
     176             :     {
     177             :       // element is local
     178        5572 :       _point_vec[0] = p;
     179        5572 :       _subproblem.reinitElemPhys(elem, _point_vec, 0);
     180        5572 :       value = _coupled_var.sln()[0];
     181             :     }
     182             :   }
     183             : 
     184             :   // broadcast value
     185        7366 :   _communicator.broadcast(value, elem_proc_id);
     186        7366 :   return value;
     187             : }
     188             : 
     189             : PostprocessorValue
     190         365 : FindValueOnLine::getValue() const
     191             : {
     192         365 :   return _position;
     193             : }