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 "FluxLimitedTVDAdvection.h"
      11             : #include "SystemBase.h"
      12             : #include "Assembly.h"
      13             : 
      14             : registerMooseObject("PorousFlowApp", FluxLimitedTVDAdvection);
      15             : 
      16             : InputParameters
      17         856 : FluxLimitedTVDAdvection::validParams()
      18             : {
      19         856 :   InputParameters params = Kernel::validParams();
      20         856 :   params.addClassDescription("Conservative form of $\\nabla \\cdot \\vec{v} u$ (advection), using "
      21             :                              "the Flux Limited TVD scheme invented by Kuzmin and Turek");
      22        1712 :   params.addRequiredParam<UserObjectName>("advective_flux_calculator",
      23             :                                           "AdvectiveFluxCalculator UserObject");
      24         856 :   return params;
      25           0 : }
      26             : 
      27         551 : FluxLimitedTVDAdvection::FluxLimitedTVDAdvection(const InputParameters & parameters)
      28             :   : Kernel(parameters),
      29         551 :     _fluo(getUserObject<AdvectiveFluxCalculatorBase>("advective_flux_calculator"))
      30             : {
      31         551 : }
      32             : 
      33             : Real
      34           0 : FluxLimitedTVDAdvection::computeQpResidual()
      35             : {
      36           0 :   mooseError("FluxLimitedTVDAdvection::computeQpResidual() called\n");
      37             :   return 0.0;
      38             : }
      39             : 
      40             : void
      41     1046348 : FluxLimitedTVDAdvection::computeResidual()
      42             : {
      43     1046348 :   prepareVectorTag(_assembly, _var.number());
      44     1046348 :   precalculateResidual();
      45             : 
      46             :   // get the residual contributions from _fluo
      47     4058484 :   for (unsigned i = 0; i < _current_elem->n_nodes(); ++i)
      48             :   {
      49     3012136 :     const dof_id_type node_id_i = _current_elem->node_id(i);
      50     3012136 :     _local_re(i) = _fluo.getFluxOut(node_id_i) / _fluo.getValence(node_id_i);
      51             :   }
      52             : 
      53     1046348 :   accumulateTaggedLocalResidual();
      54             : 
      55     1046348 :   if (_has_save_in)
      56             :   {
      57             :     Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
      58           0 :     for (unsigned int i = 0; i < _save_in.size(); i++)
      59           0 :       _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
      60             :   }
      61     1046348 : }
      62             : 
      63             : void
      64      699614 : FluxLimitedTVDAdvection::computeJacobian()
      65             : {
      66      699614 :   prepareMatrixTag(_assembly, _var.number(), _var.number());
      67      699614 :   precalculateJacobian();
      68             : 
      69             :   // Run through the nodes of this element using "i", getting the Jacobian contributions
      70             :   // d(residual_i)/du(node_j) for all nodes j that can have a nonzero Jacobian contribution.  Some
      71             :   // of these node_j will live in this element, but some will live in other elements connected with
      72             :   // node "i", and some will live in the next layer of nodes (eg, in 1D residual_3 could have
      73             :   // contributions from node1, node2, node3, node4 and node5).
      74     2658962 :   for (unsigned i = 0; i < _current_elem->n_nodes(); ++i)
      75             :   {
      76             :     // global id of node "i"
      77     1959348 :     const dof_id_type node_id_i = _current_elem->node_id(i);
      78             :     // dof number of _var on node "i"
      79             :     std::vector<dof_id_type> idof_indices(
      80     1959348 :         1, _current_elem->node_ref(i).dof_number(_sys.number(), _var.number(), 0));
      81             :     // number of times node "i" is encountered in a sweep over elements
      82     1959348 :     const unsigned valence = _fluo.getValence(node_id_i);
      83             : 
      84             :     // retrieve the derivative information from _fluo
      85     1959348 :     const std::map<dof_id_type, Real> derivs = _fluo.getdFluxOutdu(node_id_i);
      86             : 
      87             :     // now build up the dof numbers of all the "j" nodes and the derivative matrix
      88             :     // d(residual_i)/d(u_j)
      89     1959348 :     std::vector<dof_id_type> jdof_indices(derivs.size());
      90     1959348 :     DenseMatrix<Number> deriv_matrix(1, derivs.size());
      91             :     unsigned j = 0;
      92    25495680 :     for (const auto & node_j_deriv : derivs)
      93             :     {
      94             :       // global id of j:
      95    23536332 :       const dof_id_type node_id_j = node_j_deriv.first;
      96             :       // dof at node j:
      97    23536332 :       jdof_indices[j] =
      98    23536332 :           _mesh.getMesh().node_ref(node_id_j).dof_number(_sys.number(), _var.number(), 0);
      99             :       // derivative must be divided by valence, otherwise the loop over elements will multiple-count
     100    23536332 :       deriv_matrix(0, j) = node_j_deriv.second / valence;
     101    23536332 :       j++;
     102             :     }
     103             :     // Add the result to the system's Jacobian matrix
     104     1959348 :     addJacobian(_assembly, deriv_matrix, idof_indices, jdof_indices, _var.scalingFactor());
     105     3918696 :   }
     106      699614 : }