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 "PorousFlowFluxLimitedTVDAdvection.h"
      11             : #include "SystemBase.h"
      12             : #include "Assembly.h"
      13             : 
      14             : registerMooseObject("PorousFlowApp", PorousFlowFluxLimitedTVDAdvection);
      15             : 
      16             : InputParameters
      17        1796 : PorousFlowFluxLimitedTVDAdvection::validParams()
      18             : {
      19        1796 :   InputParameters params = Kernel::validParams();
      20        1796 :   params.addClassDescription("Advective flux of fluid species or heat using "
      21             :                              "the Flux Limited TVD scheme invented by Kuzmin and Turek");
      22        3592 :   params.addRequiredParam<UserObjectName>(
      23             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
      24        3592 :   params.addRequiredParam<UserObjectName>(
      25             :       "advective_flux_calculator",
      26             :       "PorousFlowAdvectiveFluxCalculator UserObject.  This determines whether the advection "
      27             :       "describes a movement of a fluid component in a fluid phase, or movement of heat energy in a "
      28             :       "fluid phase");
      29        1796 :   return params;
      30           0 : }
      31             : 
      32        1101 : PorousFlowFluxLimitedTVDAdvection::PorousFlowFluxLimitedTVDAdvection(
      33        1101 :     const InputParameters & parameters)
      34             :   : Kernel(parameters),
      35        1101 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      36        2202 :     _fluo(getUserObject<PorousFlowAdvectiveFluxCalculatorBase>("advective_flux_calculator"))
      37             : {
      38        1101 : }
      39             : 
      40             : Real
      41           0 : PorousFlowFluxLimitedTVDAdvection::computeQpResidual()
      42             : {
      43           0 :   mooseError("PorousFlowFluxLimitedTVDAdvection::computeQpResidual() called\n");
      44             :   return 0.0;
      45             : }
      46             : 
      47             : void
      48     1204142 : PorousFlowFluxLimitedTVDAdvection::computeResidual()
      49             : {
      50     1204142 :   prepareVectorTag(_assembly, _var.number());
      51     1204142 :   precalculateResidual();
      52             : 
      53             :   // get the residual contributions from _fluo
      54     3762346 :   for (unsigned i = 0; i < _current_elem->n_nodes(); ++i)
      55             :   {
      56     2558204 :     const dof_id_type node_id_i = _current_elem->node_id(i);
      57     2558204 :     _local_re(i) = _fluo.getFluxOut(node_id_i) / _fluo.getValence(node_id_i);
      58             :   }
      59             : 
      60     1204142 :   accumulateTaggedLocalResidual();
      61             : 
      62     1204142 :   if (_has_save_in)
      63             :   {
      64             :     Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
      65           0 :     for (unsigned int i = 0; i < _save_in.size(); i++)
      66           0 :       _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
      67             :   }
      68     1204142 : }
      69             : 
      70             : void
      71      891283 : PorousFlowFluxLimitedTVDAdvection::computeJacobian()
      72             : {
      73      891283 :   prepareMatrixTag(_assembly, _var.number(), _var.number());
      74      891283 :   precalculateJacobian();
      75             : 
      76             :   // Run through the nodes of this element using "i", getting the Jacobian contributions
      77             :   // d(residual_i)/du(node_j) for all nodes j that can have a nonzero Jacobian contribution.  Some
      78             :   // of these node_j will live in this element, but some will live in other elements connected with
      79             :   // node "i", and some will live in the next layer of nodes (eg, in 1D residual_3 could have
      80             :   // contributions from node1, node2, node3, node4 and node5).
      81     2769127 :   for (unsigned i = 0; i < _current_elem->n_nodes(); ++i)
      82             :   {
      83             :     // global id of node "i"
      84     1877844 :     const dof_id_type node_id_i = _current_elem->node_id(i);
      85             :     // dof number of _var on node "i"
      86             :     std::vector<dof_id_type> idof_indices(
      87     1877844 :         1, _current_elem->node_ref(i).dof_number(_sys.number(), _var.number(), 0));
      88             :     // number of times node "i" is encountered in a sweep over elements
      89     1877844 :     const unsigned valence = _fluo.getValence(node_id_i);
      90             : 
      91             :     // retrieve the derivative information from _fluo
      92     1877844 :     const std::map<dof_id_type, std::vector<Real>> derivs = _fluo.getdFluxOut_dvars(node_id_i);
      93             : 
      94             :     // now build up the dof numbers of all the "j" nodes and the derivative matrix
      95             :     // d(residual_i)/d(var_j)
      96     5644324 :     for (unsigned pvar = 0; pvar < _dictator.numVariables(); ++pvar)
      97             :     {
      98     3766480 :       const unsigned varnum = _dictator.mooseVariableNum(pvar);
      99     3766480 :       std::vector<dof_id_type> jdof_indices(derivs.size());
     100     3766480 :       DenseMatrix<Number> deriv_matrix(1, derivs.size());
     101             :       unsigned j = 0;
     102    27736936 :       for (const auto & node_j_deriv : derivs)
     103             :       {
     104             :         // global id of j:
     105    23970456 :         const dof_id_type node_id_j = node_j_deriv.first;
     106             :         // dof of pvar at node j:
     107    23970456 :         jdof_indices[j] = _mesh.getMesh().node_ref(node_id_j).dof_number(_sys.number(), varnum, 0);
     108             :         // derivative must be divided by valence, otherwise the loop over elements will
     109             :         // multiple-count
     110    23970456 :         deriv_matrix(0, j) = node_j_deriv.second[pvar] / valence;
     111    23970456 :         j++;
     112             :       }
     113             :       // Add the result to the system's Jacobian matrix
     114     3766480 :       addJacobian(_assembly, deriv_matrix, idof_indices, jdof_indices, _var.scalingFactor());
     115     3766480 :     }
     116     1877844 :   }
     117      891283 : }