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 "RichardsMultiphaseProblem.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "MooseMesh.h"
      14             : #include "MooseVariable.h"
      15             : #include "NonlinearSystem.h"
      16             : 
      17             : registerMooseObject("RichardsApp", RichardsMultiphaseProblem);
      18             : 
      19             : InputParameters
      20          10 : RichardsMultiphaseProblem::validParams()
      21             : {
      22          10 :   InputParameters params = FEProblemBase::validParams();
      23          20 :   params.addRequiredParam<NonlinearVariableName>(
      24             :       "bounded_var", "Variable whose value will be constrained to be greater than lower_var");
      25          20 :   params.addRequiredParam<NonlinearVariableName>(
      26             :       "lower_var",
      27             :       "Variable that acts as a lower bound to bounded_var.  It will not be "
      28             :       "constrained during the solution procedure");
      29          10 :   return params;
      30           0 : }
      31             : 
      32           5 : RichardsMultiphaseProblem::RichardsMultiphaseProblem(const InputParameters & params)
      33             :   : FEProblem(params),
      34             :     // in the following have to get the names of the variables, and then find their numbers in
      35             :     // initialSetup,
      36             :     // as their numbers won't be defined at the moment of instantiation of this class
      37          10 :     _bounded_var_name(params.get<NonlinearVariableName>("bounded_var")),
      38           5 :     _lower_var_name(params.get<NonlinearVariableName>("lower_var")),
      39           5 :     _bounded_var_num(0),
      40           5 :     _lower_var_num(0)
      41             : {
      42           5 : }
      43             : 
      44           6 : RichardsMultiphaseProblem::~RichardsMultiphaseProblem() {}
      45             : 
      46             : void
      47           5 : RichardsMultiphaseProblem::initialSetup()
      48             : {
      49             :   // the first argument to getVariable is threadID - i hope the following always works
      50             :   unsigned int tid = 0;
      51             : 
      52             :   // We are going to do more specific checks below, hence allowing
      53             :   // more specific error messages to be printed in case something goes
      54             :   // wrong. Therefore we just pass VAR_ANY here.
      55           5 :   MooseVariableFEBase & bounded = getVariable(
      56             :       tid, _bounded_var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD);
      57           5 :   MooseVariableFEBase & lower = getVariable(
      58             :       tid, _lower_var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD);
      59             : 
      60             :   // some checks
      61           5 :   if (!bounded.isNodal() || !lower.isNodal())
      62           1 :     mooseError("Both the bounded and lower variables must be nodal variables in "
      63             :                "RichardsMultiphaseProblem");
      64           4 :   if (bounded.feType().family != lower.feType().family)
      65           1 :     mooseError("Both the bounded and lower variables must belong to the same element family, eg "
      66             :                "LAGRANGE, in RichardsMultiphaseProblem");
      67           3 :   if (bounded.feType().order != lower.feType().order)
      68           1 :     mooseError("Both the bounded and lower variables must have the same order, eg FIRST, in "
      69             :                "RichardsMultiphaseProblem");
      70             : 
      71             :   // extract the required info
      72           2 :   _bounded_var_num = bounded.number();
      73           2 :   _lower_var_num = lower.number();
      74             : 
      75           2 :   FEProblemBase::initialSetup();
      76           2 : }
      77             : 
      78             : bool
      79          51 : RichardsMultiphaseProblem::shouldUpdateSolution()
      80             : {
      81          51 :   return true;
      82             : }
      83             : 
      84             : bool
      85          24 : RichardsMultiphaseProblem::updateSolution(NumericVector<Number> & vec_solution,
      86             :                                           NumericVector<Number> & ghosted_solution)
      87             : {
      88          24 :   bool updatedSolution =
      89             :       false; // this gets set to true if we needed to enforce the bound at any node
      90             : 
      91          24 :   unsigned int sys_num = getNonlinearSystemBase(/*nl_sys_num=*/0).number();
      92             : 
      93             :   // For parallel procs i believe that i have to use local_nodes_begin, rather than just nodes_begin
      94             :   // _mesh comes from SystemBase (_mesh = getNonlinearSystemBase().subproblem().mesh(), and
      95             :   // subproblem is this object)
      96        1056 :   for (const auto & node : _mesh.getMesh().local_node_ptr_range())
      97             :   {
      98             :     // dofs[0] is the dof number of the bounded variable at this node
      99             :     // dofs[1] is the dof number of the lower variable at this node
     100         504 :     std::vector<dof_id_type> dofs(2);
     101         504 :     dofs[0] = node->dof_number(sys_num, _bounded_var_num, 0);
     102         504 :     dofs[1] = node->dof_number(sys_num, _lower_var_num, 0);
     103             : 
     104             :     // soln[0] is the value of the bounded variable at this node
     105             :     // soln[1] is the value of the lower variable at this node
     106         504 :     std::vector<Number> soln(2);
     107         504 :     vec_solution.get(dofs, soln);
     108             : 
     109             :     // do the bounding
     110         504 :     if (soln[0] < soln[1])
     111             :     {
     112          90 :       vec_solution.set(dofs[0], soln[1]); // set the bounded variable equal to the lower value
     113          90 :       updatedSolution = true;
     114             :     }
     115         528 :   }
     116             : 
     117             :   // The above vec_solution.set calls potentially added "set" commands to a queue
     118             :   // The following actions the queue (doing MPI commands if necessary), so
     119             :   // vec_solution will actually be modified by this following command
     120          24 :   vec_solution.close();
     121             : 
     122             :   // if any proc updated the solution, all procs will know about it
     123          24 :   _communicator.max(updatedSolution);
     124             : 
     125          24 :   if (updatedSolution)
     126             :   {
     127          12 :     ghosted_solution = vec_solution;
     128          12 :     ghosted_solution.close();
     129             :   }
     130             : 
     131          24 :   return updatedSolution;
     132             : }