www.mooseframework.org
Public Member Functions | Protected Types | Protected Member Functions | Protected Attributes | List of all members
AdvectiveFluxCalculatorBase Class Referenceabstract

Base class to compute Advective fluxes. More...

#include <AdvectiveFluxCalculatorBase.h>

Inheritance diagram for AdvectiveFluxCalculatorBase:
[legend]

Public Member Functions

 AdvectiveFluxCalculatorBase (const InputParameters &parameters)
 
virtual void timestepSetup () override
 
virtual void meshChanged () override
 
virtual void initialize () override
 
virtual void threadJoin (const UserObject &uo) override
 
virtual void finalize () override
 
virtual void execute () override
 
virtual void executeOnElement (dof_id_type global_i, dof_id_type global_j, unsigned local_i, unsigned local_j, unsigned qp)
 This is called by multiple times in execute() in a double loop over _current_elem's nodes (local_i and local_j) nested in a loop over each of _current_elem's quadpoints (qp). More...
 
Real getFluxOut (dof_id_type node_i) const
 Returns the flux out of lobal node id. More...
 
const std::map< dof_id_type, Real > & getdFluxOutdu (dof_id_type node_i) const
 Returns r where r[j] = d(flux out of global node i)/du(global node j) used in Jacobian computations. More...
 
const std::vector< std::vector< Real > > & getdFluxOutdKjk (dof_id_type node_i) const
 Returns r where r[j][k] = d(flux out of global node i)/dK[connected node j][connected node k] used in Jacobian computations. More...
 
unsigned getValence (dof_id_type node_i) const
 Returns the valence of the global node i Valence is the number of times the node is encountered in a loop over elements (that have appropriate subdomain_id, if the user has employed the "blocks=" parameter) seen by this processor (including ghosted elements) More...
 

Protected Types

enum  FluxLimiterTypeEnum {
  FluxLimiterTypeEnum::MinMod, FluxLimiterTypeEnum::VanLeer, FluxLimiterTypeEnum::MC, FluxLimiterTypeEnum::superbee,
  FluxLimiterTypeEnum::None
}
 Determines Flux Limiter type (Page 135 of Kuzmin and Turek) "None" means that limitFlux=0 always, which implies zero antidiffusion will be added. More...
 
enum  PQPlusMinusEnum { PQPlusMinusEnum::PPlus, PQPlusMinusEnum::PMinus, PQPlusMinusEnum::QPlus, PQPlusMinusEnum::QMinus }
 Signals to the PQPlusMinus method what should be computed. More...
 

Protected Member Functions

virtual void buildCommLists ()
 When using multiple processors, other processors will compute: More...
 
virtual void exchangeGhostedInfo ()
 Sends and receives multi-processor information regarding u_nodal and k_ij. More...
 
virtual Real computeVelocity (unsigned i, unsigned j, unsigned qp) const =0
 Computes the transfer velocity between current node i and current node j at the current qp in the current element (_current_elem). More...
 
virtual Real computeU (unsigned i) const =0
 Computes the value of u at the local node id of the current element (_current_elem) More...
 
void limitFlux (Real a, Real b, Real &limited, Real &dlimited_db) const
 flux limiter, L, on Page 135 of Kuzmin and Turek More...
 
Real rPlus (dof_id_type sequential_i, std::vector< Real > &dlimited_du, std::vector< Real > &dlimited_dk) const
 Returns the value of R_{i}^{+}, Eqn (49) of KT. More...
 
Real rMinus (dof_id_type sequential_i, std::vector< Real > &dlimited_du, std::vector< Real > &dlimited_dk) const
 Returns the value of R_{i}^{-}, Eqn (49) of KT. More...
 
Real PQPlusMinus (dof_id_type sequential_i, const PQPlusMinusEnum pq_plus_minus, std::vector< Real > &derivs, std::vector< Real > &dpq_dk) const
 Returns the value of P_{i}^{+}, P_{i}^{-}, Q_{i}^{+} or Q_{i}^{-} (depending on pq_plus_minus) which are defined in Eqns (47) and (48) of KT. More...
 
void zeroedConnection (std::map< dof_id_type, Real > &the_map, dof_id_type node_i) const
 Clears the_map, then, using _kij, constructs the_map so that the_map[node_id] = 0.0 for all node_id connected with node_i. More...
 

Protected Attributes

bool _resizing_needed
 whether _kij, etc, need to be sized appropriately (and valence recomputed) at the start of the timestep More...
 
enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum _flux_limiter_type
 
std::vector< std::vector< Real > > _kij
 Kuzmin-Turek K_ij matrix. More...
 
std::vector< Real > _flux_out
 _flux_out[i] = flux of "heat" from sequential node i More...
 
std::vector< std::map< dof_id_type, Real > > _dflux_out_du
 _dflux_out_du[i][j] = d(flux_out[i])/d(u[j]). More...
 
std::vector< std::vector< std::vector< Real > > > _dflux_out_dKjk
 _dflux_out_dKjk[sequential_i][j][k] = d(flux_out[sequential_i])/d(K[j][k]). More...
 
std::vector< unsigned > _valence
 _valence[i] = number of times, in a loop over elements seen by this processor (viz, including ghost elements) and are part of the block-restricted blocks of this UserObject, that the sequential node i is encountered More...
 
std::vector< Real > _u_nodal
 _u_nodal[i] = value of _u at sequential node number i More...
 
std::vector< bool > _u_nodal_computed_by_thread
 _u_nodal_computed_by_thread(i) = true if _u_nodal[i] has been computed in execute() by the thread on this processor More...
 
PorousFlowConnectedNodes _connections
 Holds the sequential and global nodal IDs, and info regarding mesh connections between them. More...
 
std::size_t _number_of_nodes
 Number of nodes held by the _connections object. More...
 
processor_id_type _my_pid
 processor ID of this object More...
 
std::map< processor_id_type, std::vector< dof_id_type > > _nodes_to_receive
 _nodes_to_receive[proc_id] = list of sequential nodal IDs. More...
 
std::map< processor_id_type, std::vector< dof_id_type > > _nodes_to_send
 _nodes_to_send[proc_id] = list of sequential nodal IDs. More...
 
std::map< processor_id_type, std::vector< std::pair< dof_id_type, dof_id_type > > > _pairs_to_receive
 _pairs_to_receive[proc_id] indicates the k(i, j) pairs that will be sent to us from proc_id _pairs_to_receive is first built (in buildCommLists()) using global node IDs, but after construction, a translation to sequential node IDs and the index of connections is performed, for efficiency. More...
 
std::map< processor_id_type, std::vector< std::pair< dof_id_type, dof_id_type > > > _pairs_to_send
 _pairs_to_send[proc_id] indicates the k(i, j) pairs that we will send to proc_id _pairs_to_send is first built (in buildCommLists()) using global node IDs, but after construction, a translation to sequential node IDs and the index of connections is performed, for efficiency. More...
 
const Real _allowable_MB_wastage
 A mooseWarning is issued if mb_wasted = (_connections.sizeSequential() - _connections.numNodes()) * 4 / 1048576 > _allowable_MB_wastage. More...
 

Detailed Description

Base class to compute Advective fluxes.

Specifically, computes K_ij, D_ij, L_ij, R+, R-, f^a_ij detailed in D Kuzmin and S Turek "High-resolution FEM-TVD schemes based on a fully multidimensional flux limiter" Journal of Computational Physics 198 (2004) 131-158 Then combines the results to produce the residual and Jacobian contributions that may be used by Kernels

K_ij is a measure of flux from node i to node j D_ij is a diffusion matrix that stabilizes K_ij (K+D has the LED property) L = K + D R^+_i and R^-_i and f^a_{ij} quantify how much antidiffusion to allow around node i

Definition at line 33 of file AdvectiveFluxCalculatorBase.h.

Member Enumeration Documentation

◆ FluxLimiterTypeEnum

Determines Flux Limiter type (Page 135 of Kuzmin and Turek) "None" means that limitFlux=0 always, which implies zero antidiffusion will be added.

Enumerator
MinMod 
VanLeer 
MC 
superbee 
None 

Definition at line 169 of file AdvectiveFluxCalculatorBase.h.

169 { MinMod, VanLeer, MC, superbee, None } _flux_limiter_type;
enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum _flux_limiter_type

◆ PQPlusMinusEnum

Signals to the PQPlusMinus method what should be computed.

Enumerator
PPlus 
PMinus 
QPlus 
QMinus 

Definition at line 255 of file AdvectiveFluxCalculatorBase.h.

256  {
257  PPlus,
258  PMinus,
259  QPlus,
260  QMinus
261  };

Constructor & Destructor Documentation

◆ AdvectiveFluxCalculatorBase()

AdvectiveFluxCalculatorBase::AdvectiveFluxCalculatorBase ( const InputParameters &  parameters)

Definition at line 50 of file AdvectiveFluxCalculatorBase.C.

51  : ElementUserObject(parameters),
52  _resizing_needed(true),
53  _flux_limiter_type(getParam<MooseEnum>("flux_limiter_type").getEnum<FluxLimiterTypeEnum>()),
54  _kij(),
55  _flux_out(),
56  _dflux_out_du(),
58  _valence(),
59  _u_nodal(),
61  _connections(),
63  _my_pid(processor_id()),
68  _allowable_MB_wastage(getParam<Real>("allowable_MB_wastage"))
69 {
70  if (!_execute_enum.contains(EXEC_LINEAR))
71  paramError(
72  "execute_on",
73  "The AdvectiveFluxCalculator UserObject " + name() +
74  " execute_on parameter must include, at least, 'linear'. This is to ensure that "
75  "this UserObject computes all necessary quantities just before the Kernels evaluate "
76  "their Residuals");
77 }
std::vector< bool > _u_nodal_computed_by_thread
_u_nodal_computed_by_thread(i) = true if _u_nodal[i] has been computed in execute() by the thread on ...
std::vector< std::map< dof_id_type, Real > > _dflux_out_du
_dflux_out_du[i][j] = d(flux_out[i])/d(u[j]).
std::size_t _number_of_nodes
Number of nodes held by the _connections object.
std::vector< Real > _flux_out
_flux_out[i] = flux of "heat" from sequential node i
std::vector< std::vector< std::vector< Real > > > _dflux_out_dKjk
_dflux_out_dKjk[sequential_i][j][k] = d(flux_out[sequential_i])/d(K[j][k]).
const std::string name
Definition: Setup.h:21
enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum _flux_limiter_type
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
bool _resizing_needed
whether _kij, etc, need to be sized appropriately (and valence recomputed) at the start of the timest...
std::map< processor_id_type, std::vector< dof_id_type > > _nodes_to_send
_nodes_to_send[proc_id] = list of sequential nodal IDs.
std::map< processor_id_type, std::vector< std::pair< dof_id_type, dof_id_type > > > _pairs_to_send
_pairs_to_send[proc_id] indicates the k(i, j) pairs that we will send to proc_id _pairs_to_send is fi...
std::vector< Real > _u_nodal
_u_nodal[i] = value of _u at sequential node number i
processor_id_type _my_pid
processor ID of this object
const Real _allowable_MB_wastage
A mooseWarning is issued if mb_wasted = (_connections.sizeSequential() - _connections.numNodes()) * 4 / 1048576 > _allowable_MB_wastage.
std::vector< std::vector< Real > > _kij
Kuzmin-Turek K_ij matrix.
std::map< processor_id_type, std::vector< dof_id_type > > _nodes_to_receive
_nodes_to_receive[proc_id] = list of sequential nodal IDs.
std::map< processor_id_type, std::vector< std::pair< dof_id_type, dof_id_type > > > _pairs_to_receive
_pairs_to_receive[proc_id] indicates the k(i, j) pairs that will be sent to us from proc_id _pairs_to...
std::vector< unsigned > _valence
_valence[i] = number of times, in a loop over elements seen by this processor (viz, including ghost elements) and are part of the block-restricted blocks of this UserObject, that the sequential node i is encountered

Member Function Documentation

◆ buildCommLists()

void AdvectiveFluxCalculatorBase::buildCommLists ( )
protectedvirtual

When using multiple processors, other processors will compute:

  • u_nodal for nodes that we don't "own", but that we need when doing the stabilization
  • k_ij for node pairs that we don't "own", and contributions to node pairs that we do "own" (on the boundary of our set of elements), that are used in the stabilization This method builds _nodes_to_receive and _pairs_to_receive that describe which processors we are going to receive this info from, and similarly it builds _nodes_to_send and _pairs_to_send.

Build the multi-processor communication lists.

(A) We will have to send _u_nodal information to other processors. This is because although we can Evaluate Variables at all elements in _fe_problem.getEvaluableElementRange(), in the PorousFlow setting _u_nodal could depend on Material Properties within the elements, and we can't access those Properties within the ghosted elements.

(B) In a similar way, we need to send _kij information to other processors. A similar strategy is followed

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 868 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::buildCommLists(), and timestepSetup().

869 {
883  _nodes_to_receive.clear();
884  for (const auto & elem : _fe_problem.getEvaluableElementRange())
885  if (this->hasBlocks(elem->subdomain_id()))
886  {
887  const processor_id_type elem_pid = elem->processor_id();
888  if (elem_pid != _my_pid)
889  {
890  if (_nodes_to_receive.find(elem_pid) == _nodes_to_receive.end())
891  _nodes_to_receive[elem_pid] = std::vector<dof_id_type>();
892  for (unsigned i = 0; i < elem->n_nodes(); ++i)
893  if (std::find(_nodes_to_receive[elem_pid].begin(),
894  _nodes_to_receive[elem_pid].end(),
895  elem->node_id(i)) == _nodes_to_receive[elem_pid].end())
896  _nodes_to_receive[elem_pid].push_back(elem->node_id(i));
897  }
898  }
899 
900  // exchange this info with other processors, building _nodes_to_send at the same time
901  _nodes_to_send.clear();
902  auto nodes_action_functor = [this](processor_id_type pid, const std::vector<dof_id_type> & nts) {
903  _nodes_to_send[pid] = nts;
904  };
905  Parallel::push_parallel_vector_data(this->comm(), _nodes_to_receive, nodes_action_functor);
906 
907  // At the moment, _nodes_to_send and _nodes_to_receive contain global node numbers
908  // It is slightly more efficient to convert to sequential node numbers
909  // so that we don't have to keep doing things like
910  // _u_nodal[_connections.sequentialID(_nodes_to_send[pid][i])
911  // every time we send/receive
912  for (auto & kv : _nodes_to_send)
913  {
914  const processor_id_type pid = kv.first;
915  const std::size_t num_nodes = kv.second.size();
916  for (unsigned i = 0; i < num_nodes; ++i)
918  }
919  for (auto & kv : _nodes_to_receive)
920  {
921  const processor_id_type pid = kv.first;
922  const std::size_t num_nodes = kv.second.size();
923  for (unsigned i = 0; i < num_nodes; ++i)
925  }
926 
927  // Build pairs_to_receive
928  _pairs_to_receive.clear();
929  for (const auto & elem : _fe_problem.getEvaluableElementRange())
930  if (this->hasBlocks(elem->subdomain_id()))
931  {
932  const processor_id_type elem_pid = elem->processor_id();
933  if (elem_pid != _my_pid)
934  {
935  if (_pairs_to_receive.find(elem_pid) == _pairs_to_receive.end())
936  _pairs_to_receive[elem_pid] = std::vector<std::pair<dof_id_type, dof_id_type>>();
937  for (unsigned i = 0; i < elem->n_nodes(); ++i)
938  for (unsigned j = 0; j < elem->n_nodes(); ++j)
939  {
940  std::pair<dof_id_type, dof_id_type> the_pair(elem->node_id(i), elem->node_id(j));
941  if (std::find(_pairs_to_receive[elem_pid].begin(),
942  _pairs_to_receive[elem_pid].end(),
943  the_pair) == _pairs_to_receive[elem_pid].end())
944  _pairs_to_receive[elem_pid].push_back(the_pair);
945  }
946  }
947  }
948 
949  _pairs_to_send.clear();
950  auto pairs_action_functor = [this](processor_id_type pid,
951  const std::vector<std::pair<dof_id_type, dof_id_type>> & pts) {
952  _pairs_to_send[pid] = pts;
953  };
954  Parallel::push_parallel_vector_data(this->comm(), _pairs_to_receive, pairs_action_functor);
955 
956  // _pairs_to_send and _pairs_to_receive have been built using global node IDs
957  // since all processors know about that. However, using global IDs means
958  // that every time we send/receive, we keep having to do things like
959  // _kij[_connections.sequentialID(_pairs_to_send[pid][i].first)][_connections.indexOfGlobalConnection(_pairs_to_send[pid][i].first,
960  // _pairs_to_send[pid][i].second)] which is quite inefficient. So:
961  for (auto & kv : _pairs_to_send)
962  {
963  const processor_id_type pid = kv.first;
964  const std::size_t num_pairs = kv.second.size();
965  for (unsigned i = 0; i < num_pairs; ++i)
966  {
968  _pairs_to_send[pid][i].first, _pairs_to_send[pid][i].second);
969  _pairs_to_send[pid][i].first = _connections.sequentialID(_pairs_to_send[pid][i].first);
970  }
971  }
972  for (auto & kv : _pairs_to_receive)
973  {
974  const processor_id_type pid = kv.first;
975  const std::size_t num_pairs = kv.second.size();
976  for (unsigned i = 0; i < num_pairs; ++i)
977  {
979  _pairs_to_receive[pid][i].first, _pairs_to_receive[pid][i].second);
980  _pairs_to_receive[pid][i].first = _connections.sequentialID(_pairs_to_receive[pid][i].first);
981  }
982  }
983 }
dof_id_type sequentialID(dof_id_type global_node_ID) const
Return the sequential node ID corresponding to the global node ID This is guaranteed to lie in the ra...
unsigned indexOfGlobalConnection(dof_id_type global_node_ID_from, dof_id_type global_node_ID_to) const
Return the index of global_node_ID_to in the globalConnectionsToGlobalID(global_node_ID_from) vector...
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
std::map< processor_id_type, std::vector< dof_id_type > > _nodes_to_send
_nodes_to_send[proc_id] = list of sequential nodal IDs.
std::map< processor_id_type, std::vector< std::pair< dof_id_type, dof_id_type > > > _pairs_to_send
_pairs_to_send[proc_id] indicates the k(i, j) pairs that we will send to proc_id _pairs_to_send is fi...
processor_id_type _my_pid
processor ID of this object
std::map< processor_id_type, std::vector< dof_id_type > > _nodes_to_receive
_nodes_to_receive[proc_id] = list of sequential nodal IDs.
std::map< processor_id_type, std::vector< std::pair< dof_id_type, dof_id_type > > > _pairs_to_receive
_pairs_to_receive[proc_id] indicates the k(i, j) pairs that will be sent to us from proc_id _pairs_to...

◆ computeU()

virtual Real AdvectiveFluxCalculatorBase::computeU ( unsigned  i) const
protectedpure virtual

◆ computeVelocity()

virtual Real AdvectiveFluxCalculatorBase::computeVelocity ( unsigned  i,
unsigned  j,
unsigned  qp 
) const
protectedpure virtual

Computes the transfer velocity between current node i and current node j at the current qp in the current element (_current_elem).

For instance, (_grad_phi[i][qp] * _velocity) * _phi[j][qp];

Parameters
inode number in the current element
jnode number in the current element
qpquadpoint number in the current element

Implemented in PorousFlowAdvectiveFluxCalculatorBase, and AdvectiveFluxCalculatorConstantVelocity.

Referenced by executeOnElement().

◆ exchangeGhostedInfo()

void AdvectiveFluxCalculatorBase::exchangeGhostedInfo ( )
protectedvirtual

Sends and receives multi-processor information regarding u_nodal and k_ij.

See buildCommLists for some more explanation.

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 986 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::exchangeGhostedInfo(), and finalize().

987 {
988  // Exchange ghosted _u_nodal information with other processors
989  std::map<processor_id_type, std::vector<Real>> unodal_to_send;
990  for (const auto & kv : _nodes_to_send)
991  {
992  const processor_id_type pid = kv.first;
993  unodal_to_send[pid] = std::vector<Real>();
994  for (const auto & nd : kv.second)
995  unodal_to_send[pid].push_back(_u_nodal[nd]);
996  }
997 
998  auto unodal_action_functor = [this](processor_id_type pid,
999  const std::vector<Real> & unodal_received) {
1000  const std::size_t msg_size = unodal_received.size();
1001  mooseAssert(msg_size == _nodes_to_receive[pid].size(),
1002  "Message size, " << msg_size
1003  << ", incompatible with nodes_to_receive, which has size "
1004  << _nodes_to_receive[pid].size());
1005  for (unsigned i = 0; i < msg_size; ++i)
1006  _u_nodal[_nodes_to_receive[pid][i]] = unodal_received[i];
1007  };
1008  Parallel::push_parallel_vector_data(this->comm(), unodal_to_send, unodal_action_functor);
1009 
1010  // Exchange _kij information with other processors
1011  std::map<processor_id_type, std::vector<Real>> kij_to_send;
1012  for (const auto & kv : _pairs_to_send)
1013  {
1014  const processor_id_type pid = kv.first;
1015  kij_to_send[pid] = std::vector<Real>();
1016  for (const auto & pr : kv.second)
1017  kij_to_send[pid].push_back(_kij[pr.first][pr.second]);
1018  }
1019 
1020  auto kij_action_functor = [this](processor_id_type pid, const std::vector<Real> & kij_received) {
1021  const std::size_t msg_size = kij_received.size();
1022  mooseAssert(msg_size == _pairs_to_receive[pid].size(),
1023  "Message size, " << msg_size
1024  << ", incompatible with pairs_to_receive, which has size "
1025  << _pairs_to_receive[pid].size());
1026  for (unsigned i = 0; i < msg_size; ++i)
1027  _kij[_pairs_to_receive[pid][i].first][_pairs_to_receive[pid][i].second] += kij_received[i];
1028  };
1029  Parallel::push_parallel_vector_data(this->comm(), kij_to_send, kij_action_functor);
1030 }
std::map< processor_id_type, std::vector< dof_id_type > > _nodes_to_send
_nodes_to_send[proc_id] = list of sequential nodal IDs.
std::map< processor_id_type, std::vector< std::pair< dof_id_type, dof_id_type > > > _pairs_to_send
_pairs_to_send[proc_id] indicates the k(i, j) pairs that we will send to proc_id _pairs_to_send is fi...
std::vector< Real > _u_nodal
_u_nodal[i] = value of _u at sequential node number i
std::vector< std::vector< Real > > _kij
Kuzmin-Turek K_ij matrix.
std::map< processor_id_type, std::vector< dof_id_type > > _nodes_to_receive
_nodes_to_receive[proc_id] = list of sequential nodal IDs.
std::map< processor_id_type, std::vector< std::pair< dof_id_type, dof_id_type > > > _pairs_to_receive
_pairs_to_receive[proc_id] indicates the k(i, j) pairs that will be sent to us from proc_id _pairs_to...

◆ execute()

void AdvectiveFluxCalculatorBase::execute ( )
overridevirtual

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 199 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::execute().

200 {
201  // compute _kij contributions from this element that is local to this processor
202  // and record _u_nodal
203  for (unsigned i = 0; i < _current_elem->n_nodes(); ++i)
204  {
205  const dof_id_type node_i = _current_elem->node_id(i);
206  const dof_id_type sequential_i = _connections.sequentialID(node_i);
207  if (!_u_nodal_computed_by_thread[sequential_i])
208  {
209  _u_nodal[sequential_i] = computeU(i);
210  _u_nodal_computed_by_thread[sequential_i] = true;
211  }
212  for (unsigned j = 0; j < _current_elem->n_nodes(); ++j)
213  {
214  const dof_id_type node_j = _current_elem->node_id(j);
215  for (unsigned qp = 0; qp < _qrule->n_points(); ++qp)
216  executeOnElement(node_i, node_j, i, j, qp);
217  }
218  }
219 }
std::vector< bool > _u_nodal_computed_by_thread
_u_nodal_computed_by_thread(i) = true if _u_nodal[i] has been computed in execute() by the thread on ...
dof_id_type sequentialID(dof_id_type global_node_ID) const
Return the sequential node ID corresponding to the global node ID This is guaranteed to lie in the ra...
virtual void executeOnElement(dof_id_type global_i, dof_id_type global_j, unsigned local_i, unsigned local_j, unsigned qp)
This is called by multiple times in execute() in a double loop over _current_elem&#39;s nodes (local_i an...
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
std::vector< Real > _u_nodal
_u_nodal[i] = value of _u at sequential node number i
virtual Real computeU(unsigned i) const =0
Computes the value of u at the local node id of the current element (_current_elem) ...

◆ executeOnElement()

void AdvectiveFluxCalculatorBase::executeOnElement ( dof_id_type  global_i,
dof_id_type  global_j,
unsigned  local_i,
unsigned  local_j,
unsigned  qp 
)
virtual

This is called by multiple times in execute() in a double loop over _current_elem's nodes (local_i and local_j) nested in a loop over each of _current_elem's quadpoints (qp).

It is used to compute _kij and its derivatives

Parameters
global_iglobal node id corresponding to the local node local_i
global_jglobal node id corresponding to the local node local_j
local_ilocal node number of the _current_elem
local_jlocal node number of the _current_elem
qpquadpoint number of the _current_elem

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 222 of file AdvectiveFluxCalculatorBase.C.

Referenced by execute(), and PorousFlowAdvectiveFluxCalculatorBase::executeOnElement().

224 {
225  // KT Eqn (20)
226  const dof_id_type sequential_i = _connections.sequentialID(global_i);
227  const unsigned index_i_to_j = _connections.indexOfGlobalConnection(global_i, global_j);
228  _kij[sequential_i][index_i_to_j] += _JxW[qp] * _coord[qp] * computeVelocity(local_i, local_j, qp);
229 }
dof_id_type sequentialID(dof_id_type global_node_ID) const
Return the sequential node ID corresponding to the global node ID This is guaranteed to lie in the ra...
unsigned indexOfGlobalConnection(dof_id_type global_node_ID_from, dof_id_type global_node_ID_to) const
Return the index of global_node_ID_to in the globalConnectionsToGlobalID(global_node_ID_from) vector...
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
virtual Real computeVelocity(unsigned i, unsigned j, unsigned qp) const =0
Computes the transfer velocity between current node i and current node j at the current qp in the cur...
std::vector< std::vector< Real > > _kij
Kuzmin-Turek K_ij matrix.

◆ finalize()

void AdvectiveFluxCalculatorBase::finalize ( )
overridevirtual

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 251 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::finalize().

252 {
253  // Overall: ensure _kij is fully built, then compute Kuzmin-Turek D, L, f^a and
254  // relevant Jacobian information, and then the relevant quantities into _flux_out and
255  // _dflux_out_du, _dflux_out_dKjk
256 
257  if (_app.n_processors() > 1)
259 
260  // Calculate KuzminTurek D matrix
261  // See Eqn (32)
262  // This adds artificial diffusion, which eliminates any spurious oscillations
263  // The idea is that D will remove all negative off-diagonal elements when it is added to K
264  // This is identical to full upwinding
265  std::vector<std::vector<Real>> dij(_number_of_nodes);
266  std::vector<std::vector<Real>> dDij_dKij(
267  _number_of_nodes); // dDij_dKij[i][j] = d(D[i][j])/d(K[i][j]) for i!=j
268  std::vector<std::vector<Real>> dDij_dKji(
269  _number_of_nodes); // dDij_dKji[i][j] = d(D[i][j])/d(K[j][i]) for i!=j
270  std::vector<std::vector<Real>> dDii_dKij(
271  _number_of_nodes); // dDii_dKij[i][j] = d(D[i][i])/d(K[i][j])
272  std::vector<std::vector<Real>> dDii_dKji(
273  _number_of_nodes); // dDii_dKji[i][j] = d(D[i][i])/d(K[j][i])
274  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
275  {
276  const std::vector<dof_id_type> con_i =
278  const std::size_t num_con_i = con_i.size();
279  dij[sequential_i].assign(num_con_i, 0.0);
280  dDij_dKij[sequential_i].assign(num_con_i, 0.0);
281  dDij_dKji[sequential_i].assign(num_con_i, 0.0);
282  dDii_dKij[sequential_i].assign(num_con_i, 0.0);
283  dDii_dKji[sequential_i].assign(num_con_i, 0.0);
284  const unsigned index_i_to_i =
285  _connections.indexOfSequentialConnection(sequential_i, sequential_i);
286  for (std::size_t j = 0; j < num_con_i; ++j)
287  {
288  const dof_id_type sequential_j = con_i[j];
289  if (sequential_i == sequential_j)
290  continue;
291  const unsigned index_j_to_i =
292  _connections.indexOfSequentialConnection(sequential_j, sequential_i);
293  const Real kij = _kij[sequential_i][j];
294  const Real kji = _kij[sequential_j][index_j_to_i];
295  if ((kij <= kji) && (kij < 0.0))
296  {
297  dij[sequential_i][j] = -kij;
298  dDij_dKij[sequential_i][j] = -1.0;
299  dDii_dKij[sequential_i][j] += 1.0;
300  }
301  else if ((kji <= kij) && (kji < 0.0))
302  {
303  dij[sequential_i][j] = -kji;
304  dDij_dKji[sequential_i][j] = -1.0;
305  dDii_dKji[sequential_i][j] += 1.0;
306  }
307  else
308  dij[sequential_i][j] = 0.0;
309  dij[sequential_i][index_i_to_i] -= dij[sequential_i][j];
310  }
311  }
312 
313  // Calculate KuzminTurek L matrix
314  // See Fig 2: L = K + D
315  std::vector<std::vector<Real>> lij(_number_of_nodes);
316  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
317  {
318  const std::vector<dof_id_type> con_i =
320  const std::size_t num_con_i = con_i.size();
321  lij[sequential_i].assign(num_con_i, 0.0);
322  for (std::size_t j = 0; j < num_con_i; ++j)
323  lij[sequential_i][j] = _kij[sequential_i][j] + dij[sequential_i][j];
324  }
325 
326  // Compute KuzminTurek R matrices
327  // See Eqns (49) and (12)
328  std::vector<Real> rP(_number_of_nodes);
329  std::vector<Real> rM(_number_of_nodes);
330  std::vector<std::vector<Real>> drP(_number_of_nodes); // drP[i][j] = d(rP[i])/d(u[j]). Here j
331  // indexes the j^th node connected to i
332  std::vector<std::vector<Real>> drM(_number_of_nodes);
333  std::vector<std::vector<Real>> drP_dk(
334  _number_of_nodes); // drP_dk[i][j] = d(rP[i])/d(K[i][j]). Here j indexes the j^th node
335  // connected to i
336  std::vector<std::vector<Real>> drM_dk(_number_of_nodes);
337  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
338  {
339  rP[sequential_i] = rPlus(sequential_i, drP[sequential_i], drP_dk[sequential_i]);
340  rM[sequential_i] = rMinus(sequential_i, drM[sequential_i], drM_dk[sequential_i]);
341  }
342 
343  // Calculate KuzminTurek f^{a} matrix
344  // This is the antidiffusive flux
345  // See Eqn (50)
346  std::vector<std::vector<Real>> fa(_number_of_nodes); // fa[sequential_i][j] sequential_j is the
347  // j^th connection to sequential_i
348  // The derivatives are a bit complicated.
349  // If i is upwind of j then fa[i][j] depends on all nodes connected to i.
350  // But if i is downwind of j then fa[i][j] depends on all nodes connected to j.
351  std::vector<std::vector<std::map<dof_id_type, Real>>> dfa(
352  _number_of_nodes); // dfa[sequential_i][j][global_k] = d(fa[sequential_i][j])/du[global_k].
353  // Here global_k can be a neighbor to sequential_i or a neighbour to
354  // sequential_j (sequential_j is the j^th connection to sequential_i)
355  std::vector<std::vector<std::vector<Real>>> dFij_dKik(
356  _number_of_nodes); // dFij_dKik[sequential_i][j][k] =
357  // d(fa[sequential_i][j])/d(K[sequential_i][k]). Here j denotes the j^th
358  // connection to sequential_i, while k denotes the k^th connection to
359  // sequential_i
360  std::vector<std::vector<std::vector<Real>>> dFij_dKjk(
361  _number_of_nodes); // dFij_dKjk[sequential_i][j][k] =
362  // d(fa[sequential_i][j])/d(K[sequential_j][k]). Here sequential_j is
363  // the j^th connection to sequential_i, and k denotes the k^th connection
364  // to sequential_j (this will include sequential_i itself)
365  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
366  {
367  const std::vector<dof_id_type> con_i =
369  const unsigned num_con_i = con_i.size();
370  fa[sequential_i].assign(num_con_i, 0.0);
371  dfa[sequential_i].resize(num_con_i);
372  dFij_dKik[sequential_i].resize(num_con_i);
373  dFij_dKjk[sequential_i].resize(num_con_i);
374  for (std::size_t j = 0; j < num_con_i; ++j)
375  {
376  for (const auto & global_k : con_i)
377  dfa[sequential_i][j][global_k] = 0;
378  const dof_id_type global_j = con_i[j];
379  const std::vector<dof_id_type> con_j = _connections.globalConnectionsToGlobalID(global_j);
380  const unsigned num_con_j = con_j.size();
381  for (const auto & global_k : con_j)
382  dfa[sequential_i][j][global_k] = 0;
383  dFij_dKik[sequential_i][j].assign(num_con_i, 0.0);
384  dFij_dKjk[sequential_i][j].assign(num_con_j, 0.0);
385  }
386  }
387 
388  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
389  {
390  const dof_id_type global_i = _connections.globalID(sequential_i);
391  const Real u_i = _u_nodal[sequential_i];
392  const std::vector<dof_id_type> con_i =
394  const std::size_t num_con_i = con_i.size();
395  for (std::size_t j = 0; j < num_con_i; ++j)
396  {
397  const dof_id_type sequential_j = con_i[j];
398  if (sequential_i == sequential_j)
399  continue;
400  const unsigned index_j_to_i =
401  _connections.indexOfSequentialConnection(sequential_j, sequential_i);
402  const Real Lij = lij[sequential_i][j];
403  const Real Lji = lij[sequential_j][index_j_to_i];
404  if (Lji >= Lij) // node i is upwind of node j.
405  {
406  const Real Dij = dij[sequential_i][j];
407  const dof_id_type global_j = _connections.globalID(sequential_j);
408  const Real u_j = _u_nodal[sequential_j];
409  Real prefactor = 0.0;
410  std::vector<Real> dprefactor_du(num_con_i,
411  0.0); // dprefactor_du[j] = d(prefactor)/du[sequential_j]
412  std::vector<Real> dprefactor_dKij(
413  num_con_i, 0.0); // dprefactor_dKij[j] = d(prefactor)/dKij[sequential_i][j].
414  Real dprefactor_dKji = 0; // dprefactor_dKji = d(prefactor)/dKij[sequential_j][index_j_to_i]
415  if (u_i >= u_j)
416  {
417  if (Lji <= rP[sequential_i] * Dij)
418  {
419  prefactor = Lji;
420  dprefactor_dKij[j] += dDij_dKji[sequential_j][index_j_to_i];
421  dprefactor_dKji += 1.0 + dDij_dKij[sequential_j][index_j_to_i];
422  }
423  else
424  {
425  prefactor = rP[sequential_i] * Dij;
426  for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
427  dprefactor_du[ind_j] = drP[sequential_i][ind_j] * Dij;
428  dprefactor_dKij[j] += rP[sequential_i] * dDij_dKij[sequential_i][j];
429  dprefactor_dKji += rP[sequential_i] * dDij_dKji[sequential_i][j];
430  for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
431  dprefactor_dKij[ind_j] += drP_dk[sequential_i][ind_j] * Dij;
432  }
433  }
434  else
435  {
436  if (Lji <= rM[sequential_i] * Dij)
437  {
438  prefactor = Lji;
439  dprefactor_dKij[j] += dDij_dKji[sequential_j][index_j_to_i];
440  dprefactor_dKji += 1.0 + dDij_dKij[sequential_j][index_j_to_i];
441  }
442  else
443  {
444  prefactor = rM[sequential_i] * Dij;
445  for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
446  dprefactor_du[ind_j] = drM[sequential_i][ind_j] * Dij;
447  dprefactor_dKij[j] += rM[sequential_i] * dDij_dKij[sequential_i][j];
448  dprefactor_dKji += rM[sequential_i] * dDij_dKji[sequential_i][j];
449  for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
450  dprefactor_dKij[ind_j] += drM_dk[sequential_i][ind_j] * Dij;
451  }
452  }
453  fa[sequential_i][j] = prefactor * (u_i - u_j);
454  dfa[sequential_i][j][global_i] = prefactor;
455  dfa[sequential_i][j][global_j] = -prefactor;
456  for (std::size_t ind_j = 0; ind_j < num_con_i; ++ind_j)
457  {
458  dfa[sequential_i][j][_connections.globalID(con_i[ind_j])] +=
459  dprefactor_du[ind_j] * (u_i - u_j);
460  dFij_dKik[sequential_i][j][ind_j] += dprefactor_dKij[ind_j] * (u_i - u_j);
461  }
462  dFij_dKjk[sequential_i][j][index_j_to_i] += dprefactor_dKji * (u_i - u_j);
463  }
464  }
465  }
466 
467  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
468  {
469  const std::vector<dof_id_type> con_i =
471  const std::size_t num_con_i = con_i.size();
472  for (std::size_t j = 0; j < num_con_i; ++j)
473  {
474  const dof_id_type sequential_j = con_i[j];
475  if (sequential_i == sequential_j)
476  continue;
477  const unsigned index_j_to_i =
478  _connections.indexOfSequentialConnection(sequential_j, sequential_i);
479  if (lij[sequential_j][index_j_to_i] < lij[sequential_i][j]) // node_i is downwind of node_j.
480  {
481  fa[sequential_i][j] = -fa[sequential_j][index_j_to_i];
482  for (const auto & dof_deriv : dfa[sequential_j][index_j_to_i])
483  dfa[sequential_i][j][dof_deriv.first] = -dof_deriv.second;
484  for (std::size_t k = 0; k < num_con_i; ++k)
485  dFij_dKik[sequential_i][j][k] = -dFij_dKjk[sequential_j][index_j_to_i][k];
486  const std::size_t num_con_j =
488  for (std::size_t k = 0; k < num_con_j; ++k)
489  dFij_dKjk[sequential_i][j][k] = -dFij_dKik[sequential_j][index_j_to_i][k];
490  }
491  }
492  }
493 
494  // zero _flux_out and its derivatives
495  _flux_out.assign(_number_of_nodes, 0.0);
496  // The derivatives are a bit complicated.
497  // If i is upwind of a node "j" then _flux_out[i] depends on all nodes connected to i.
498  // But if i is downwind of a node "j" then _flux_out depends on all nodes connected with node
499  // j.
500  _dflux_out_du.assign(_number_of_nodes, std::map<dof_id_type, Real>());
501  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
502  for (const auto & j : _connections.globalConnectionsToSequentialID(sequential_i))
503  {
504  _dflux_out_du[sequential_i][j] = 0.0;
505  for (const auto & neighbors_j : _connections.globalConnectionsToGlobalID(j))
506  _dflux_out_du[sequential_i][neighbors_j] = 0.0;
507  }
509  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
510  {
511  const std::vector<dof_id_type> con_i =
513  const std::size_t num_con_i = con_i.size();
514  _dflux_out_dKjk[sequential_i].resize(num_con_i);
515  for (std::size_t j = 0; j < num_con_i; ++j)
516  {
517  const dof_id_type sequential_j = con_i[j];
518  std::vector<dof_id_type> con_j =
520  _dflux_out_dKjk[sequential_i][j].assign(con_j.size(), 0.0);
521  }
522  }
523 
524  // Add everything together
525  // See step 3 in Fig 2, noting Eqn (36)
526  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
527  {
528  std::vector<dof_id_type> con_i = _connections.sequentialConnectionsToSequentialID(sequential_i);
529  const size_t num_con_i = con_i.size();
530  const dof_id_type index_i_to_i =
531  _connections.indexOfSequentialConnection(sequential_i, sequential_i);
532  for (std::size_t j = 0; j < num_con_i; ++j)
533  {
534  const dof_id_type sequential_j = con_i[j];
535  const dof_id_type global_j = _connections.globalID(sequential_j);
536  const Real u_j = _u_nodal[sequential_j];
537 
538  // negative sign because residual = -Lu (KT equation (19))
539  _flux_out[sequential_i] -= lij[sequential_i][j] * u_j + fa[sequential_i][j];
540 
541  _dflux_out_du[sequential_i][global_j] -= lij[sequential_i][j];
542  for (const auto & dof_deriv : dfa[sequential_i][j])
543  _dflux_out_du[sequential_i][dof_deriv.first] -= dof_deriv.second;
544 
545  _dflux_out_dKjk[sequential_i][index_i_to_i][j] -= 1.0 * u_j; // from the K in L = K + D
546 
547  if (sequential_j == sequential_i)
548  for (dof_id_type k = 0; k < num_con_i; ++k)
549  _dflux_out_dKjk[sequential_i][index_i_to_i][k] -= dDii_dKij[sequential_i][k] * u_j;
550  else
551  _dflux_out_dKjk[sequential_i][index_i_to_i][j] -= dDij_dKij[sequential_i][j] * u_j;
552  for (dof_id_type k = 0; k < num_con_i; ++k)
553  _dflux_out_dKjk[sequential_i][index_i_to_i][k] -= dFij_dKik[sequential_i][j][k];
554 
555  if (sequential_j == sequential_i)
556  for (unsigned k = 0; k < con_i.size(); ++k)
557  {
558  const unsigned index_k_to_i =
559  _connections.indexOfSequentialConnection(con_i[k], sequential_i);
560  _dflux_out_dKjk[sequential_i][k][index_k_to_i] -= dDii_dKji[sequential_i][k] * u_j;
561  }
562  else
563  {
564  const unsigned index_j_to_i =
565  _connections.indexOfSequentialConnection(sequential_j, sequential_i);
566  _dflux_out_dKjk[sequential_i][j][index_j_to_i] -= dDij_dKji[sequential_i][j] * u_j;
567  }
568  for (unsigned k = 0;
569  k < _connections.sequentialConnectionsToSequentialID(sequential_j).size();
570  ++k)
571  _dflux_out_dKjk[sequential_i][j][k] -= dFij_dKjk[sequential_i][j][k];
572  }
573  }
574 }
dof_id_type globalID(dof_id_type sequential_node_ID) const
Return the global node ID (node number in the mesh) corresponding to the provided sequential node ID...
std::vector< std::map< dof_id_type, Real > > _dflux_out_du
_dflux_out_du[i][j] = d(flux_out[i])/d(u[j]).
const std::vector< dof_id_type > & sequentialConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (sequential node IDs) connected to the given sequential node ID All elements of ...
Real rMinus(dof_id_type sequential_i, std::vector< Real > &dlimited_du, std::vector< Real > &dlimited_dk) const
Returns the value of R_{i}^{-}, Eqn (49) of KT.
std::size_t _number_of_nodes
Number of nodes held by the _connections object.
std::vector< Real > _flux_out
_flux_out[i] = flux of "heat" from sequential node i
std::vector< std::vector< std::vector< Real > > > _dflux_out_dKjk
_dflux_out_dKjk[sequential_i][j][k] = d(flux_out[sequential_i])/d(K[j][k]).
const std::vector< dof_id_type > & globalConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (global node IDs) connected to the given sequential node ID. ...
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
Real rPlus(dof_id_type sequential_i, std::vector< Real > &dlimited_du, std::vector< Real > &dlimited_dk) const
Returns the value of R_{i}^{+}, Eqn (49) of KT.
std::vector< Real > _u_nodal
_u_nodal[i] = value of _u at sequential node number i
virtual void exchangeGhostedInfo()
Sends and receives multi-processor information regarding u_nodal and k_ij.
std::vector< std::vector< Real > > _kij
Kuzmin-Turek K_ij matrix.
const std::vector< dof_id_type > & globalConnectionsToGlobalID(dof_id_type global_node_ID) const
Return all the nodes (global node IDs) connected to the given global node ID.
unsigned indexOfSequentialConnection(dof_id_type sequential_node_ID_from, dof_id_type sequential_node_ID_to) const
Return the index of sequential_node_ID_to in the sequentialConnectionsToSequentialID(sequential_node_...

◆ getdFluxOutdKjk()

const std::vector< std::vector< Real > > & AdvectiveFluxCalculatorBase::getdFluxOutdKjk ( dof_id_type  node_i) const

Returns r where r[j][k] = d(flux out of global node i)/dK[connected node j][connected node k] used in Jacobian computations.

Parameters
node_iglobal id of node
Returns
the derivatives (after applying the KT procedure)

Definition at line 755 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::finalize().

756 {
757  return _dflux_out_dKjk[_connections.sequentialID(node_i)];
758 }
dof_id_type sequentialID(dof_id_type global_node_ID) const
Return the sequential node ID corresponding to the global node ID This is guaranteed to lie in the ra...
std::vector< std::vector< std::vector< Real > > > _dflux_out_dKjk
_dflux_out_dKjk[sequential_i][j][k] = d(flux_out[sequential_i])/d(K[j][k]).
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...

◆ getdFluxOutdu()

const std::map< dof_id_type, Real > & AdvectiveFluxCalculatorBase::getdFluxOutdu ( dof_id_type  node_i) const

Returns r where r[j] = d(flux out of global node i)/du(global node j) used in Jacobian computations.

Parameters
node_iglobal id of node
Returns
the derivatives (after applying the KT procedure)

Definition at line 749 of file AdvectiveFluxCalculatorBase.C.

Referenced by FluxLimitedTVDAdvection::computeJacobian(), and PorousFlowAdvectiveFluxCalculatorBase::finalize().

750 {
751  return _dflux_out_du[_connections.sequentialID(node_i)];
752 }
std::vector< std::map< dof_id_type, Real > > _dflux_out_du
_dflux_out_du[i][j] = d(flux_out[i])/d(u[j]).
dof_id_type sequentialID(dof_id_type global_node_ID) const
Return the sequential node ID corresponding to the global node ID This is guaranteed to lie in the ra...
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...

◆ getFluxOut()

Real AdvectiveFluxCalculatorBase::getFluxOut ( dof_id_type  node_i) const

Returns the flux out of lobal node id.

Parameters
node_iid of node
Returns
advective flux out of node after applying the KT procedure

Definition at line 761 of file AdvectiveFluxCalculatorBase.C.

Referenced by FluxLimitedTVDAdvection::computeResidual(), and PorousFlowFluxLimitedTVDAdvection::computeResidual().

762 {
763  return _flux_out[_connections.sequentialID(node_i)];
764 }
dof_id_type sequentialID(dof_id_type global_node_ID) const
Return the sequential node ID corresponding to the global node ID This is guaranteed to lie in the ra...
std::vector< Real > _flux_out
_flux_out[i] = flux of "heat" from sequential node i
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...

◆ getValence()

unsigned AdvectiveFluxCalculatorBase::getValence ( dof_id_type  node_i) const

Returns the valence of the global node i Valence is the number of times the node is encountered in a loop over elements (that have appropriate subdomain_id, if the user has employed the "blocks=" parameter) seen by this processor (including ghosted elements)

Parameters
node_igloal id of i^th node
Returns
valence of the node

Definition at line 767 of file AdvectiveFluxCalculatorBase.C.

Referenced by FluxLimitedTVDAdvection::computeJacobian(), PorousFlowFluxLimitedTVDAdvection::computeJacobian(), FluxLimitedTVDAdvection::computeResidual(), and PorousFlowFluxLimitedTVDAdvection::computeResidual().

768 {
769  return _valence[_connections.sequentialID(node_i)];
770 }
dof_id_type sequentialID(dof_id_type global_node_ID) const
Return the sequential node ID corresponding to the global node ID This is guaranteed to lie in the ra...
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
std::vector< unsigned > _valence
_valence[i] = number of times, in a loop over elements seen by this processor (viz, including ghost elements) and are part of the block-restricted blocks of this UserObject, that the sequential node i is encountered

◆ initialize()

void AdvectiveFluxCalculatorBase::initialize ( )
overridevirtual

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 188 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::initialize().

189 {
190  // Zero _kij and falsify _u_nodal_computed_by_thread ready for building in execute() and
191  // finalize()
193  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
194  _kij[sequential_i].assign(_connections.sequentialConnectionsToSequentialID(sequential_i).size(),
195  0.0);
196 }
std::vector< bool > _u_nodal_computed_by_thread
_u_nodal_computed_by_thread(i) = true if _u_nodal[i] has been computed in execute() by the thread on ...
const std::vector< dof_id_type > & sequentialConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (sequential node IDs) connected to the given sequential node ID All elements of ...
std::size_t _number_of_nodes
Number of nodes held by the _connections object.
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
std::vector< std::vector< Real > > _kij
Kuzmin-Turek K_ij matrix.

◆ limitFlux()

void AdvectiveFluxCalculatorBase::limitFlux ( Real  a,
Real  b,
Real &  limited,
Real &  dlimited_db 
) const
protected

flux limiter, L, on Page 135 of Kuzmin and Turek

Parameters
aKT's "a" parameter
bKT's "b" parameter
limited[out]The value of the flux limiter, L
dlimited_db[out]The derivative dL/db

Definition at line 649 of file AdvectiveFluxCalculatorBase.C.

Referenced by rMinus(), and rPlus().

650 {
651  limited = 0.0;
652  dlimited_db = 0.0;
654  return;
655 
656  if ((a >= 0.0 && b <= 0.0) || (a <= 0.0 && b >= 0.0))
657  return;
658  const Real s = (a > 0.0 ? 1.0 : -1.0);
659 
660  const Real lal = std::abs(a);
661  const Real lbl = std::abs(b);
662  const Real dlbl = (b >= 0.0 ? 1.0 : -1.0); // d(lbl)/db
663  switch (_flux_limiter_type)
664  {
666  {
667  if (lal <= lbl)
668  {
669  limited = s * lal;
670  dlimited_db = 0.0;
671  }
672  else
673  {
674  limited = s * lbl;
675  dlimited_db = s * dlbl;
676  }
677  return;
678  }
680  {
681  limited = s * 2 * lal * lbl / (lal + lbl);
682  dlimited_db = s * 2 * lal * (dlbl / (lal + lbl) - lbl * dlbl / std::pow(lal + lbl, 2));
683  return;
684  }
686  {
687  const Real av = 0.5 * std::abs(a + b);
688  if (2 * lal <= av && lal <= lbl)
689  {
690  // 2 * lal is the smallest
691  limited = s * 2.0 * lal;
692  dlimited_db = 0.0;
693  }
694  else if (2 * lbl <= av && lbl <= lal)
695  {
696  // 2 * lbl is the smallest
697  limited = s * 2.0 * lbl;
698  dlimited_db = s * 2.0 * dlbl;
699  }
700  else
701  {
702  // av is the smallest
703  limited = s * av;
704  // if (a>0 and b>0) then d(av)/db = 0.5 = 0.5 * dlbl
705  // if (a<0 and b<0) then d(av)/db = -0.5 = 0.5 * dlbl
706  // if a and b have different sign then limited=0, above
707  dlimited_db = s * 0.5 * dlbl;
708  }
709  return;
710  }
712  {
713  const Real term1 = std::min(2.0 * lal, lbl);
714  const Real term2 = std::min(lal, 2.0 * lbl);
715  if (term1 >= term2)
716  {
717  if (2.0 * lal <= lbl)
718  {
719  limited = s * 2 * lal;
720  dlimited_db = 0.0;
721  }
722  else
723  {
724  limited = s * lbl;
725  dlimited_db = s * dlbl;
726  }
727  }
728  else
729  {
730  if (lal <= 2.0 * lbl)
731  {
732  limited = s * lal;
733  dlimited_db = 0.0;
734  }
735  else
736  {
737  limited = s * 2.0 * lbl;
738  dlimited_db = s * 2.0 * dlbl;
739  }
740  }
741  return;
742  }
743  default:
744  return;
745  }
746 }
enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum _flux_limiter_type
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)

◆ meshChanged()

void AdvectiveFluxCalculatorBase::meshChanged ( )
overridevirtual

Definition at line 179 of file AdvectiveFluxCalculatorBase.C.

180 {
181  ElementUserObject::meshChanged();
182 
183  // Signal that _kij, _valence, etc need to be rebuilt
184  _resizing_needed = true;
185 }
bool _resizing_needed
whether _kij, etc, need to be sized appropriately (and valence recomputed) at the start of the timest...

◆ PQPlusMinus()

Real AdvectiveFluxCalculatorBase::PQPlusMinus ( dof_id_type  sequential_i,
const PQPlusMinusEnum  pq_plus_minus,
std::vector< Real > &  derivs,
std::vector< Real > &  dpq_dk 
) const
protected

Returns the value of P_{i}^{+}, P_{i}^{-}, Q_{i}^{+} or Q_{i}^{-} (depending on pq_plus_minus) which are defined in Eqns (47) and (48) of KT.

Parameters
sequential_isequential nodal ID
pq_plus_minusindicates whether P_{i}^{+}, P_{i}^{-}, Q_{i}^{+} or Q_{i}^{-} should be returned
derivs[out]derivs[j] = d(result)/d(u[sequential_j]). Here sequential_j is the j^th connection to sequential_i
dpq_dk[out]dpq_dk[j] = d(result)/d(K[node_i][j]). Here j indexes a connection to sequential_i. Recall that d(result)/d(K[l][m]) are zero unless l=sequential_i

Definition at line 782 of file AdvectiveFluxCalculatorBase.C.

Referenced by rMinus(), and rPlus().

786 {
787  // Find the value of u at sequential_i
788  const Real u_i = _u_nodal[sequential_i];
789 
790  // Connections to sequential_i
791  const std::vector<dof_id_type> con_i =
793  const std::size_t num_con = con_i.size();
794  // The neighbor number of sequential_i to sequential_i
795  const unsigned i_index_i = _connections.indexOfSequentialConnection(sequential_i, sequential_i);
796 
797  // Initialize the results
798  Real result = 0.0;
799  derivs.assign(num_con, 0.0);
800  dpqdk.assign(num_con, 0.0);
801 
802  // Sum over all nodes connected with node_i.
803  for (std::size_t j = 0; j < num_con; ++j)
804  {
805  const dof_id_type sequential_j = con_i[j];
806  if (sequential_j == sequential_i)
807  continue;
808  const Real kentry = _kij[sequential_i][j];
809 
810  // Find the value of u at node_j
811  const Real u_j = _u_nodal[sequential_j];
812  const Real ujminusi = u_j - u_i;
813 
814  // Evaluate the i-j contribution to the result
815  switch (pq_plus_minus)
816  {
818  {
819  if (ujminusi < 0.0 && kentry < 0.0)
820  {
821  result += kentry * ujminusi;
822  derivs[j] += kentry;
823  derivs[i_index_i] -= kentry;
824  dpqdk[j] += ujminusi;
825  }
826  break;
827  }
829  {
830  if (ujminusi > 0.0 && kentry < 0.0)
831  {
832  result += kentry * ujminusi;
833  derivs[j] += kentry;
834  derivs[i_index_i] -= kentry;
835  dpqdk[j] += ujminusi;
836  }
837  break;
838  }
840  {
841  if (ujminusi > 0.0 && kentry > 0.0)
842  {
843  result += kentry * ujminusi;
844  derivs[j] += kentry;
845  derivs[i_index_i] -= kentry;
846  dpqdk[j] += ujminusi;
847  }
848  break;
849  }
851  {
852  if (ujminusi < 0.0 && kentry > 0.0)
853  {
854  result += kentry * ujminusi;
855  derivs[j] += kentry;
856  derivs[i_index_i] -= kentry;
857  dpqdk[j] += ujminusi;
858  }
859  break;
860  }
861  }
862  }
863 
864  return result;
865 }
const std::vector< dof_id_type > & sequentialConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (sequential node IDs) connected to the given sequential node ID All elements of ...
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
std::vector< Real > _u_nodal
_u_nodal[i] = value of _u at sequential node number i
std::vector< std::vector< Real > > _kij
Kuzmin-Turek K_ij matrix.
unsigned indexOfSequentialConnection(dof_id_type sequential_node_ID_from, dof_id_type sequential_node_ID_to) const
Return the index of sequential_node_ID_to in the sequentialConnectionsToSequentialID(sequential_node_...

◆ rMinus()

Real AdvectiveFluxCalculatorBase::rMinus ( dof_id_type  sequential_i,
std::vector< Real > &  dlimited_du,
std::vector< Real > &  dlimited_dk 
) const
protected

Returns the value of R_{i}^{-}, Eqn (49) of KT.

Parameters
sequential_iSequential nodal ID
dlimited_du[out]dlimited_du[j] = d(R_{sequential_i}^{-})/du[sequential_j]. Here sequential_j is the j^th connection to sequential_i
dlimited_dk[out]dlimited_dk[j] = d(R_{sequential_i}^{-})/d(K[sequential_i][j]). Note Derivatives w.r.t. K[l][m] with l!=i are zero

Definition at line 613 of file AdvectiveFluxCalculatorBase.C.

Referenced by finalize().

616 {
617  const std::size_t num_con = _connections.sequentialConnectionsToSequentialID(sequential_i).size();
618  dlimited_du.assign(num_con, 0.0);
619  dlimited_dk.assign(num_con, 0.0);
621  return 0.0;
622  std::vector<Real> dp_du;
623  std::vector<Real> dp_dk;
624  const Real p = PQPlusMinus(sequential_i, PQPlusMinusEnum::PMinus, dp_du, dp_dk);
625  if (p == 0.0)
626  // Comment after Eqn (49): if P=0 then there's no antidiffusion, so no need to remove it
627  return 1.0;
628  std::vector<Real> dq_du;
629  std::vector<Real> dq_dk;
630  const Real q = PQPlusMinus(sequential_i, PQPlusMinusEnum::QMinus, dq_du, dq_dk);
631 
632  const Real r = q / p;
633  Real limited;
634  Real dlimited_dr;
635  limitFlux(1.0, r, limited, dlimited_dr);
636 
637  const Real p2 = std::pow(p, 2);
638  for (std::size_t j = 0; j < num_con; ++j)
639  {
640  const Real dr_du = dq_du[j] / p - q * dp_du[j] / p2;
641  const Real dr_dk = dq_dk[j] / p - q * dp_dk[j] / p2;
642  dlimited_du[j] = dlimited_dr * dr_du;
643  dlimited_dk[j] = dlimited_dr * dr_dk;
644  }
645  return limited;
646 }
void limitFlux(Real a, Real b, Real &limited, Real &dlimited_db) const
flux limiter, L, on Page 135 of Kuzmin and Turek
const std::vector< dof_id_type > & sequentialConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (sequential node IDs) connected to the given sequential node ID All elements of ...
Real PQPlusMinus(dof_id_type sequential_i, const PQPlusMinusEnum pq_plus_minus, std::vector< Real > &derivs, std::vector< Real > &dpq_dk) const
Returns the value of P_{i}^{+}, P_{i}^{-}, Q_{i}^{+} or Q_{i}^{-} (depending on pq_plus_minus) which ...
enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum _flux_limiter_type
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)

◆ rPlus()

Real AdvectiveFluxCalculatorBase::rPlus ( dof_id_type  sequential_i,
std::vector< Real > &  dlimited_du,
std::vector< Real > &  dlimited_dk 
) const
protected

Returns the value of R_{i}^{+}, Eqn (49) of KT.

Parameters
node_inodal id
dlimited_du[out]dlimited_du[j] = d(R_{sequential_i}^{+})/du[sequential_j]. Here sequential_j is the j^th connection to sequential_i
dlimited_dk[out]dlimited_dk[j] = d(R_{sequential_i}^{+})/d(K[sequential_i][j]). Note Derivatives w.r.t. K[l][m] with l!=i are zero

Definition at line 577 of file AdvectiveFluxCalculatorBase.C.

Referenced by finalize().

580 {
581  const std::size_t num_con = _connections.sequentialConnectionsToSequentialID(sequential_i).size();
582  dlimited_du.assign(num_con, 0.0);
583  dlimited_dk.assign(num_con, 0.0);
585  return 0.0;
586  std::vector<Real> dp_du;
587  std::vector<Real> dp_dk;
588  const Real p = PQPlusMinus(sequential_i, PQPlusMinusEnum::PPlus, dp_du, dp_dk);
589  if (p == 0.0)
590  // Comment after Eqn (49): if P=0 then there's no antidiffusion, so no need to remove it
591  return 1.0;
592  std::vector<Real> dq_du;
593  std::vector<Real> dq_dk;
594  const Real q = PQPlusMinus(sequential_i, PQPlusMinusEnum::QPlus, dq_du, dq_dk);
595 
596  const Real r = q / p;
597  Real limited;
598  Real dlimited_dr;
599  limitFlux(1.0, r, limited, dlimited_dr);
600 
601  const Real p2 = std::pow(p, 2);
602  for (std::size_t j = 0; j < num_con; ++j)
603  {
604  const Real dr_du = dq_du[j] / p - q * dp_du[j] / p2;
605  const Real dr_dk = dq_dk[j] / p - q * dp_dk[j] / p2;
606  dlimited_du[j] = dlimited_dr * dr_du;
607  dlimited_dk[j] = dlimited_dr * dr_dk;
608  }
609  return limited;
610 }
void limitFlux(Real a, Real b, Real &limited, Real &dlimited_db) const
flux limiter, L, on Page 135 of Kuzmin and Turek
const std::vector< dof_id_type > & sequentialConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (sequential node IDs) connected to the given sequential node ID All elements of ...
Real PQPlusMinus(dof_id_type sequential_i, const PQPlusMinusEnum pq_plus_minus, std::vector< Real > &derivs, std::vector< Real > &dpq_dk) const
Returns the value of P_{i}^{+}, P_{i}^{-}, Q_{i}^{+} or Q_{i}^{-} (depending on pq_plus_minus) which ...
enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum _flux_limiter_type
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)

◆ threadJoin()

void AdvectiveFluxCalculatorBase::threadJoin ( const UserObject &  uo)
overridevirtual

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 232 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::threadJoin().

233 {
234  const AdvectiveFluxCalculatorBase & afc = static_cast<const AdvectiveFluxCalculatorBase &>(uo);
235  // add the values of _kij computed by different threads
236  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
237  {
238  const std::size_t num_con_i =
240  for (std::size_t j = 0; j < num_con_i; ++j)
241  _kij[sequential_i][j] += afc._kij[sequential_i][j];
242  }
243 
244  // gather the values of _u_nodal computed by different threads
245  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
246  if (!_u_nodal_computed_by_thread[sequential_i] && afc._u_nodal_computed_by_thread[sequential_i])
247  _u_nodal[sequential_i] = afc._u_nodal[sequential_i];
248 }
std::vector< bool > _u_nodal_computed_by_thread
_u_nodal_computed_by_thread(i) = true if _u_nodal[i] has been computed in execute() by the thread on ...
const std::vector< dof_id_type > & sequentialConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (sequential node IDs) connected to the given sequential node ID All elements of ...
std::size_t _number_of_nodes
Number of nodes held by the _connections object.
Base class to compute Advective fluxes.
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
std::vector< Real > _u_nodal
_u_nodal[i] = value of _u at sequential node number i
std::vector< std::vector< Real > > _kij
Kuzmin-Turek K_ij matrix.

◆ timestepSetup()

void AdvectiveFluxCalculatorBase::timestepSetup ( )
overridevirtual

Reimplemented in PorousFlowAdvectiveFluxCalculatorBase.

Definition at line 80 of file AdvectiveFluxCalculatorBase.C.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::timestepSetup().

81 {
82  // If needed, size and initialize quantities appropriately, and compute _valence
83  if (_resizing_needed)
84  {
85  /*
86  * Populate _connections for all nodes that can be seen by this processor and on relevant
87  * blocks
88  *
89  * MULTIPROC NOTE: this must loop over local elements and 2 layers of ghosted elements.
90  * The Kernel will only loop over local elements, so will only use _kij, etc, for
91  * linked node-node pairs that appear in the local elements. Nevertheless, we
92  * need to build _kij, etc, for the nodes in the ghosted elements in order to simplify
93  * Jacobian computations
94  */
96  for (const auto & elem : _fe_problem.getEvaluableElementRange())
97  if (this->hasBlocks(elem->subdomain_id()))
98  for (unsigned i = 0; i < elem->n_nodes(); ++i)
99  _connections.addGlobalNode(elem->node_id(i));
101  for (const auto & elem : _fe_problem.getEvaluableElementRange())
102  if (this->hasBlocks(elem->subdomain_id()))
103  for (unsigned i = 0; i < elem->n_nodes(); ++i)
104  for (unsigned j = 0; j < elem->n_nodes(); ++j)
105  _connections.addConnection(elem->node_id(i), elem->node_id(j));
107 
109 
110  Real mb_wasted = (_connections.sizeSequential() - _number_of_nodes) * 4.0 / 1048576;
111  gatherMax(mb_wasted);
112  if (mb_wasted > _allowable_MB_wastage)
113  mooseWarning(
114  "In at least one processor, the sequential node-numbering internal data structure used "
115  "in " +
116  name() + " is using memory inefficiently.\nThe memory wasted is " +
117  Moose::stringify(mb_wasted) +
118  " megabytes.\n The node-numbering data structure has been optimized for speed at the "
119  "expense of memory, but that may not be an appropriate optimization for your case, "
120  "because the node numbering is not close to sequential in your case.\n");
121 
122  // initialize _kij
123  _kij.resize(_number_of_nodes);
124  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
125  _kij[sequential_i].assign(
126  _connections.sequentialConnectionsToSequentialID(sequential_i).size(), 0.0);
127 
128  /*
129  * Build _valence[i], which is the number of times the sequential node i is encountered when
130  * looping over the entire mesh (and on relevant blocks)
131  *
132  * MULTIPROC NOTE: this must loop over local elements and >=1 layer of ghosted elements.
133  * The Kernel will only loop over local elements, so will only use _valence for
134  * linked node-node pairs that appear in the local elements. But other processors will
135  * loop over neighboring elements, so avoid multiple counting of the residual and Jacobian
136  * this processor must record how many times each node-node link of its local elements
137  * appears in the local elements and >=1 layer, and pass that info to the Kernel
138  */
139  _valence.assign(_number_of_nodes, 0);
140  for (const auto & elem : _fe_problem.getEvaluableElementRange())
141  if (this->hasBlocks(elem->subdomain_id()))
142  for (unsigned i = 0; i < elem->n_nodes(); ++i)
143  {
144  const dof_id_type node_i = elem->node_id(i);
145  const dof_id_type sequential_i = _connections.sequentialID(node_i);
146  _valence[sequential_i] += 1;
147  }
148 
149  _u_nodal.assign(_number_of_nodes, 0.0);
151  _flux_out.assign(_number_of_nodes, 0.0);
152  _dflux_out_du.assign(_number_of_nodes, std::map<dof_id_type, Real>());
153  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
154  zeroedConnection(_dflux_out_du[sequential_i], _connections.globalID(sequential_i));
156  for (dof_id_type sequential_i = 0; sequential_i < _number_of_nodes; ++sequential_i)
157  {
158  const std::vector<dof_id_type> con_i =
160  const std::size_t num_con_i = con_i.size();
161  _dflux_out_dKjk[sequential_i].resize(num_con_i);
162  for (std::size_t j = 0; j < con_i.size(); ++j)
163  {
164  const dof_id_type sequential_j = con_i[j];
165  const std::size_t num_con_j =
167  _dflux_out_dKjk[sequential_i][j].assign(num_con_j, 0.0);
168  }
169  }
170 
171  if (_app.n_processors() > 1)
172  buildCommLists();
173 
174  _resizing_needed = false;
175  }
176 }
void clear()
clear all data in readiness for adding global nodes and connections
std::vector< bool > _u_nodal_computed_by_thread
_u_nodal_computed_by_thread(i) = true if _u_nodal[i] has been computed in execute() by the thread on ...
dof_id_type globalID(dof_id_type sequential_node_ID) const
Return the global node ID (node number in the mesh) corresponding to the provided sequential node ID...
std::vector< std::map< dof_id_type, Real > > _dflux_out_du
_dflux_out_du[i][j] = d(flux_out[i])/d(u[j]).
void addConnection(dof_id_type global_node_from, dof_id_type global_node_to)
Specifies that global_node_to is connected to global_node_from.
dof_id_type sequentialID(dof_id_type global_node_ID) const
Return the sequential node ID corresponding to the global node ID This is guaranteed to lie in the ra...
const std::vector< dof_id_type > & sequentialConnectionsToSequentialID(dof_id_type sequential_node_ID) const
Return all the nodes (sequential node IDs) connected to the given sequential node ID All elements of ...
void finalizeAddingConnections()
Signal that all global node IDs have been added to the internal data structures.
std::size_t _number_of_nodes
Number of nodes held by the _connections object.
std::vector< Real > _flux_out
_flux_out[i] = flux of "heat" from sequential node i
std::vector< std::vector< std::vector< Real > > > _dflux_out_dKjk
_dflux_out_dKjk[sequential_i][j][k] = d(flux_out[sequential_i])/d(K[j][k]).
std::size_t sizeSequential() const
Return the size of _sequential_id, for checking memory efficiency.
const std::string name
Definition: Setup.h:21
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
bool _resizing_needed
whether _kij, etc, need to be sized appropriately (and valence recomputed) at the start of the timest...
std::vector< Real > _u_nodal
_u_nodal[i] = value of _u at sequential node number i
const Real _allowable_MB_wastage
A mooseWarning is issued if mb_wasted = (_connections.sizeSequential() - _connections.numNodes()) * 4 / 1048576 > _allowable_MB_wastage.
std::size_t numNodes() const
number of nodes known by this class
virtual void buildCommLists()
When using multiple processors, other processors will compute:
std::vector< std::vector< Real > > _kij
Kuzmin-Turek K_ij matrix.
void finalizeAddingGlobalNodes()
Signal that all global node IDs have been added to the internal data structures.
void zeroedConnection(std::map< dof_id_type, Real > &the_map, dof_id_type node_i) const
Clears the_map, then, using _kij, constructs the_map so that the_map[node_id] = 0.0 for all node_id connected with node_i.
std::vector< unsigned > _valence
_valence[i] = number of times, in a loop over elements seen by this processor (viz, including ghost elements) and are part of the block-restricted blocks of this UserObject, that the sequential node i is encountered
void addGlobalNode(dof_id_type global_node_ID)
Add the given global_node_ID to the internal data structures If the global node ID has already been a...

◆ zeroedConnection()

void AdvectiveFluxCalculatorBase::zeroedConnection ( std::map< dof_id_type, Real > &  the_map,
dof_id_type  node_i 
) const
protected

Clears the_map, then, using _kij, constructs the_map so that the_map[node_id] = 0.0 for all node_id connected with node_i.

Parameters
[out]the_mapthe map to be zeroed appropriately
[in]node_inodal id

Definition at line 773 of file AdvectiveFluxCalculatorBase.C.

Referenced by timestepSetup().

775 {
776  the_map.clear();
777  for (const auto & node_j : _connections.globalConnectionsToGlobalID(node_i))
778  the_map[node_j] = 0.0;
779 }
PorousFlowConnectedNodes _connections
Holds the sequential and global nodal IDs, and info regarding mesh connections between them...
const std::vector< dof_id_type > & globalConnectionsToGlobalID(dof_id_type global_node_ID) const
Return all the nodes (global node IDs) connected to the given global node ID.

Member Data Documentation

◆ _allowable_MB_wastage

const Real AdvectiveFluxCalculatorBase::_allowable_MB_wastage
protected

A mooseWarning is issued if mb_wasted = (_connections.sizeSequential() - _connections.numNodes()) * 4 / 1048576 > _allowable_MB_wastage.

The _connections object uses sequential node numbering for computational efficiency, but this leads to memory being used inefficiently: the number of megabytes wasted is mb_wasted.

Definition at line 252 of file AdvectiveFluxCalculatorBase.h.

Referenced by timestepSetup().

◆ _connections

PorousFlowConnectedNodes AdvectiveFluxCalculatorBase::_connections
protected

◆ _dflux_out_dKjk

std::vector<std::vector<std::vector<Real> > > AdvectiveFluxCalculatorBase::_dflux_out_dKjk
protected

_dflux_out_dKjk[sequential_i][j][k] = d(flux_out[sequential_i])/d(K[j][k]).

Here sequential_i is a sequential node number according to the _connections object, and j represents the j^th node connected to i according to the _connections object, and k represents the k^th node connected to j according to the _connections object. Here j must be connected to i (this does include (the sequential version of j) == i), and k must be connected to j (this does include (the sequential version of k) = i and (the sequential version of k) == (sequential version of j)))

Definition at line 190 of file AdvectiveFluxCalculatorBase.h.

Referenced by finalize(), getdFluxOutdKjk(), and timestepSetup().

◆ _dflux_out_du

std::vector<std::map<dof_id_type, Real> > AdvectiveFluxCalculatorBase::_dflux_out_du
protected

_dflux_out_du[i][j] = d(flux_out[i])/d(u[j]).

Here i is a sequential node number according to the _connections object, and j (global ID) must be connected to i, or to a node that is connected to i.

Definition at line 183 of file AdvectiveFluxCalculatorBase.h.

Referenced by finalize(), getdFluxOutdu(), and timestepSetup().

◆ _flux_limiter_type

enum AdvectiveFluxCalculatorBase::FluxLimiterTypeEnum AdvectiveFluxCalculatorBase::_flux_limiter_type
protected

Referenced by limitFlux(), rMinus(), and rPlus().

◆ _flux_out

std::vector<Real> AdvectiveFluxCalculatorBase::_flux_out
protected

_flux_out[i] = flux of "heat" from sequential node i

Definition at line 179 of file AdvectiveFluxCalculatorBase.h.

Referenced by finalize(), getFluxOut(), and timestepSetup().

◆ _kij

std::vector<std::vector<Real> > AdvectiveFluxCalculatorBase::_kij
protected

Kuzmin-Turek K_ij matrix.

Along with R+ and R-, this is the key quantity computed by this UserObject. _kij[i][j] = k_ij corresponding to the i-j node pair. Here i is a sequential node numbers according to the _connections object, and j represents the j^th node connected to i according to the _connections object.

Definition at line 176 of file AdvectiveFluxCalculatorBase.h.

Referenced by exchangeGhostedInfo(), executeOnElement(), finalize(), initialize(), PQPlusMinus(), threadJoin(), and timestepSetup().

◆ _my_pid

processor_id_type AdvectiveFluxCalculatorBase::_my_pid
protected

processor ID of this object

Definition at line 210 of file AdvectiveFluxCalculatorBase.h.

Referenced by PorousFlowAdvectiveFluxCalculatorBase::buildCommLists(), and buildCommLists().

◆ _nodes_to_receive

std::map<processor_id_type, std::vector<dof_id_type> > AdvectiveFluxCalculatorBase::_nodes_to_receive
protected

_nodes_to_receive[proc_id] = list of sequential nodal IDs.

proc_id will send us _u_nodal at those nodes. _nodes_to_receive is built (in buildCommLists()) using global node IDs, but after construction, a translation to sequential node IDs is made, for efficiency. The result is: we will receive _u_nodal[_nodes_to_receive[proc_id][i]] from proc_id

Definition at line 218 of file AdvectiveFluxCalculatorBase.h.

Referenced by buildCommLists(), PorousFlowAdvectiveFluxCalculatorBase::exchangeGhostedInfo(), and exchangeGhostedInfo().

◆ _nodes_to_send

std::map<processor_id_type, std::vector<dof_id_type> > AdvectiveFluxCalculatorBase::_nodes_to_send
protected

_nodes_to_send[proc_id] = list of sequential nodal IDs.

We will send _u_nodal at those nodes to proc_id _nodes_to_send is built (in buildCommLists()) using global node IDs, but after construction, a translation to sequential node IDs is made, for efficiency The result is: we will send _u_nodal[_nodes_to_receive[proc_id][i]] to proc_id

Definition at line 226 of file AdvectiveFluxCalculatorBase.h.

Referenced by buildCommLists(), PorousFlowAdvectiveFluxCalculatorBase::exchangeGhostedInfo(), and exchangeGhostedInfo().

◆ _number_of_nodes

std::size_t AdvectiveFluxCalculatorBase::_number_of_nodes
protected

Number of nodes held by the _connections object.

Definition at line 207 of file AdvectiveFluxCalculatorBase.h.

Referenced by finalize(), initialize(), threadJoin(), PorousFlowAdvectiveFluxCalculatorBase::threadJoin(), and timestepSetup().

◆ _pairs_to_receive

std::map<processor_id_type, std::vector<std::pair<dof_id_type, dof_id_type> > > AdvectiveFluxCalculatorBase::_pairs_to_receive
protected

_pairs_to_receive[proc_id] indicates the k(i, j) pairs that will be sent to us from proc_id _pairs_to_receive is first built (in buildCommLists()) using global node IDs, but after construction, a translation to sequential node IDs and the index of connections is performed, for efficiency.

The result is we will receive: _kij[_pairs_to_receive[proc_id][i].first][_pairs_to_receive[proc_id][i].second] from proc_id

Definition at line 235 of file AdvectiveFluxCalculatorBase.h.

Referenced by buildCommLists(), and exchangeGhostedInfo().

◆ _pairs_to_send

std::map<processor_id_type, std::vector<std::pair<dof_id_type, dof_id_type> > > AdvectiveFluxCalculatorBase::_pairs_to_send
protected

_pairs_to_send[proc_id] indicates the k(i, j) pairs that we will send to proc_id _pairs_to_send is first built (in buildCommLists()) using global node IDs, but after construction, a translation to sequential node IDs and the index of connections is performed, for efficiency.

The result is we will send: _kij[_pairs_to_send[proc_id][i].first][_pairs_to_send[proc_id][i+1].second] to proc_id

Definition at line 244 of file AdvectiveFluxCalculatorBase.h.

Referenced by buildCommLists(), and exchangeGhostedInfo().

◆ _resizing_needed

bool AdvectiveFluxCalculatorBase::_resizing_needed
protected

whether _kij, etc, need to be sized appropriately (and valence recomputed) at the start of the timestep

Definition at line 130 of file AdvectiveFluxCalculatorBase.h.

Referenced by meshChanged(), timestepSetup(), and PorousFlowAdvectiveFluxCalculatorBase::timestepSetup().

◆ _u_nodal

std::vector<Real> AdvectiveFluxCalculatorBase::_u_nodal
protected

_u_nodal[i] = value of _u at sequential node number i

Definition at line 198 of file AdvectiveFluxCalculatorBase.h.

Referenced by exchangeGhostedInfo(), execute(), finalize(), PQPlusMinus(), threadJoin(), and timestepSetup().

◆ _u_nodal_computed_by_thread

std::vector<bool> AdvectiveFluxCalculatorBase::_u_nodal_computed_by_thread
protected

_u_nodal_computed_by_thread(i) = true if _u_nodal[i] has been computed in execute() by the thread on this processor

Definition at line 201 of file AdvectiveFluxCalculatorBase.h.

Referenced by execute(), initialize(), threadJoin(), and timestepSetup().

◆ _valence

std::vector<unsigned> AdvectiveFluxCalculatorBase::_valence
protected

_valence[i] = number of times, in a loop over elements seen by this processor (viz, including ghost elements) and are part of the block-restricted blocks of this UserObject, that the sequential node i is encountered

Definition at line 195 of file AdvectiveFluxCalculatorBase.h.

Referenced by getValence(), and timestepSetup().


The documentation for this class was generated from the following files: