NonlinearThread.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#include "NonlinearThread.h"
#include "AuxiliarySystem.h"
#include "NonlinearSystem.h"
#include "Problem.h"
#include "FEProblem.h"
#include "KernelBase.h"
#include "IntegratedBCBase.h"
#include "DGKernelBase.h"
#include "InterfaceKernelBase.h"
#include "Material.h"
#include "TimeKernel.h"
#include "SwapBackSentinel.h"
#include "FVTimeKernel.h"
#include "ComputeJacobianThread.h"
#include "Assembly.h"
 
#include "libmesh/threads.h"
#include "libmesh/elem.h"
#include "libmesh/fe_map.h"
 
NonlinearThread::NonlinearThread(FEProblemBase & fe_problem)
  : ThreadedElementLoop<ConstElemRange>(fe_problem),
    _nl(fe_problem.currentNonlinearSystem()),
    _num_cached(0),
    _integrated_bcs(_nl.getIntegratedBCWarehouse()),
    _dg_kernels(_nl.getDGKernelWarehouse()),
    _interface_kernels(_nl.getInterfaceKernelWarehouse()),
    _kernels(_nl.getKernelWarehouse()),
    _hdg_kernels(_nl.getHDGKernelWarehouse()),
    _has_active_objects(_integrated_bcs.hasActiveObjects() || _dg_kernels.hasActiveObjects() ||
                        _interface_kernels.hasActiveObjects() || _kernels.hasActiveObjects() ||
                        _fe_problem.haveFV()),
    _should_execute_dg(false)
{
}
 
// Splitting Constructor
NonlinearThread::NonlinearThread(NonlinearThread & x, Threads::split split)
  : ThreadedElementLoop<ConstElemRange>(x, split),
    _nl(x._nl),
    _num_cached(x._num_cached),
    _integrated_bcs(x._integrated_bcs),
    _dg_kernels(x._dg_kernels),
    _interface_kernels(x._interface_kernels),
    _kernels(x._kernels),
    _tag_kernels(x._tag_kernels),
    _hdg_kernels(x._hdg_kernels),
    _has_active_objects(x._has_active_objects),
    _should_execute_dg(x._should_execute_dg)
{
}
 
NonlinearThread::~NonlinearThread() {}
 
void
NonlinearThread::operator()(const ConstElemRange & range, bool bypass_threading)
{
  if (_has_active_objects)
    ThreadedElementLoop<ConstElemRange>::operator()(range, bypass_threading);
}
 
void
NonlinearThread::subdomainChanged()
{
  // This should come first to setup the residual objects before we do dependency determination of
  // material properties and variables
  determineObjectWarehouses();
 
  _fe_problem.subdomainSetup(_subdomain, _tid);
 
  // Update variable Dependencies
  std::set<MooseVariableFEBase *> needed_moose_vars;
  _tag_kernels->updateBlockVariableDependency(_subdomain, needed_moose_vars, _tid);
  _ibc_warehouse->updateBoundaryVariableDependency(needed_moose_vars, _tid);
  _dg_warehouse->updateBlockVariableDependency(_subdomain, needed_moose_vars, _tid);
  _ik_warehouse->updateBoundaryVariableDependency(needed_moose_vars, _tid);
 
  // Update FE variable coupleable vector tags
  std::set<TagID> needed_fe_var_vector_tags;
  _tag_kernels->updateBlockFEVariableCoupledVectorTagDependency(
      _subdomain, needed_fe_var_vector_tags, _tid);
  _ibc_warehouse->updateBlockFEVariableCoupledVectorTagDependency(
      _subdomain, needed_fe_var_vector_tags, _tid);
  _fe_problem.getMaterialWarehouse().updateBlockFEVariableCoupledVectorTagDependency(
      _subdomain, needed_fe_var_vector_tags, _tid);
 
  // Update material dependencies
  std::unordered_set<unsigned int> needed_mat_props;
  _tag_kernels->updateBlockMatPropDependency(_subdomain, needed_mat_props, _tid);
  _ibc_warehouse->updateBoundaryMatPropDependency(needed_mat_props, _tid);
  _dg_warehouse->updateBlockMatPropDependency(_subdomain, needed_mat_props, _tid);
  _ik_warehouse->updateBoundaryMatPropDependency(needed_mat_props, _tid);
 
  if (_fe_problem.haveFV())
  {
    // Re-query the finite volume elemental kernels
    _fv_kernels.clear();
    _fe_problem.theWarehouse()
        .query()
        .template condition<AttribSysNum>(_nl.number())
        .template condition<AttribSystem>("FVElementalKernel")
        .template condition<AttribSubdomains>(_subdomain)
        .template condition<AttribThread>(_tid)
        .queryInto(_fv_kernels);
    for (const auto fv_kernel : _fv_kernels)
    {
      const auto & fv_mv_deps = fv_kernel->getMooseVariableDependencies();
      needed_moose_vars.insert(fv_mv_deps.begin(), fv_mv_deps.end());
      const auto & fv_mp_deps = fv_kernel->getMatPropDependencies();
      needed_mat_props.insert(fv_mp_deps.begin(), fv_mp_deps.end());
    }
  }
 
  // Cache these to avoid computing them on every side
  _subdomain_has_dg = _dg_warehouse->hasActiveBlockObjects(_subdomain, _tid);
  _subdomain_has_hdg = _hdg_warehouse->hasActiveBlockObjects(_subdomain, _tid);
 
  _fe_problem.setActiveElementalMooseVariables(needed_moose_vars, _tid);
  _fe_problem.setActiveFEVariableCoupleableVectorTags(needed_fe_var_vector_tags, _tid);
  _fe_problem.prepareMaterials(needed_mat_props, _subdomain, _tid);
}
 
void
NonlinearThread::onElement(const Elem * const elem)
{
  // Set up Sentinel class so that, even if reinitMaterials() throws in prepareElement, we
  // still remember to swap back during stack unwinding.
  SwapBackSentinel sentinel(_fe_problem, &FEProblemBase::swapBackMaterials, this->_tid);
 
  prepareElement(elem);
 
  if (dynamic_cast<ComputeJacobianThread *>(this))
    if (_nl.getScalarVariables(_tid).size() > 0)
      _fe_problem.reinitOffDiagScalars(_tid);
 
  computeOnElement();
}
 
void
NonlinearThread::computeOnElement()
{
  if (_tag_kernels->hasActiveBlockObjects(_subdomain, _tid))
  {
    const auto & kernels = _tag_kernels->getActiveBlockObjects(_subdomain, _tid);
    for (const auto & kernel : kernels)
      compute(*kernel);
  }
 
  if (_fe_problem.haveFV())
    for (auto kernel : _fv_kernels)
      compute(*kernel);
}
 
void
NonlinearThread::onBoundary(const Elem * const elem,
                            const unsigned int side,
                            const BoundaryID bnd_id,
                            const Elem * const lower_d_elem /*=nullptr*/)
{
  if (_ibc_warehouse->hasActiveBoundaryObjects(bnd_id, _tid))
  {
    // Set up Sentinel class so that, after we swap in reinitMaterialsFace in prepareFace, even if
    // one of our callees throws we remember to swap back during stack unwinding. We put our
    // sentinel here as opposed to in prepareFace because we certainly don't want our materials
    // swapped back before we proceed to residual/Jacobian computation
    SwapBackSentinel sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
 
    prepareFace(elem, side, bnd_id, lower_d_elem);
    computeOnBoundary(bnd_id, lower_d_elem);
 
    if (lower_d_elem)
      accumulateLower();
  }
}
 
void
NonlinearThread::computeOnBoundary(BoundaryID bnd_id, const Elem * /*lower_d_elem*/)
{
  const auto & bcs = _ibc_warehouse->getActiveBoundaryObjects(bnd_id, _tid);
  for (const auto & bc : bcs)
    if (bc->shouldApply())
      compute(*bc);
}
 
void
NonlinearThread::onInterface(const Elem * elem, unsigned int side, BoundaryID bnd_id)
{
  if (_ik_warehouse->hasActiveBoundaryObjects(bnd_id, _tid))
  {
 
    // Pointer to the neighbor we are currently working on.
    const Elem * neighbor = elem->neighbor_ptr(side);
 
    if (neighbor->active())
    {
      _fe_problem.reinitNeighbor(elem, side, _tid);
 
      // Set up Sentinels so that, even if one of the reinitMaterialsXXX() calls throws, we
      // still remember to swap back during stack unwinding. Note that face, boundary, and interface
      // all operate with the same MaterialData object
      SwapBackSentinel face_sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
      _fe_problem.reinitMaterialsFaceOnBoundary(bnd_id, elem->subdomain_id(), _tid);
      _fe_problem.reinitMaterialsBoundary(bnd_id, _tid);
 
      SwapBackSentinel neighbor_sentinel(_fe_problem, &FEProblem::swapBackMaterialsNeighbor, _tid);
      _fe_problem.reinitMaterialsNeighborOnBoundary(bnd_id, neighbor->subdomain_id(), _tid);
 
      // Has to happen after face and neighbor properties have been computed. Note that we don't use
      // a sentinel here because FEProblem::swapBackMaterialsFace is going to handle face materials,
      // boundary materials, and interface materials (e.g. it queries the boundary material data
      // with the current element and side
      _fe_problem.reinitMaterialsInterface(bnd_id, _tid);
 
      computeOnInterface(bnd_id);
 
      accumulateNeighbor();
    }
    else
    {
      // elem's neighbor across this side is an inactive ancestor: elem is on the coarser side of
      // an interface with h-refinement on the other side.
      // The real geometric neighbors are the active children of this ancestor, so visit each of
      // them individually instead of silently skipping the whole face.
      std::vector<const Elem *> fine_neighbors;
      neighbor->active_family_tree_by_neighbor(fine_neighbors, elem);
      for (const Elem * fine_neighbor : fine_neighbors)
        onInterfaceWithHRefinementMismatch(elem, side, bnd_id, fine_neighbor);
    }
  }
}
 
void
NonlinearThread::onInterfaceWithHRefinementMismatch(const Elem * elem,
                                                    unsigned int side,
                                                    BoundaryID bnd_id,
                                                    const Elem * fine_neighbor)
{
  if (_fe_problem.getDisplacedProblem())
    mooseError("InterfaceKernel '",
               _ik_warehouse->getActiveBoundaryObjects(bnd_id, _tid)[0]->name(),
               "' is defined on boundary ",
               bnd_id,
               ", whose primary (boundary-tagged) side is coarser than its paired side. This "
               "case is not supported with a displaced mesh problem: the displaced problem's "
               "assembly is not reinitialized on these faces, and the borrowed integration "
               "weight used to handle the mismatched refinement is not AD-differentiable with "
               "respect to mesh displacement.");
 
  const unsigned int fine_side = fine_neighbor->which_neighbor_am_i(elem);
 
  // Measure the fine neighbor's own natural face quadrature: these physical points and weights
  // are the correct integration measure for the portion of elem's face that fine_neighbor
  // actually covers. Call this directly on the Assembly object (not through FEProblemBase) so no
  // variable/dof binding happens for fine_neighbor in the "elem" role.
  Assembly & current_assembly = _fe_problem.assembly(_tid, _fe_problem.currentNlSysNum());
  current_assembly.reinitElemFaceRef(fine_neighbor, fine_side, TOLERANCE, nullptr, nullptr);
  const auto fine_phys_points = current_assembly.qPointsFace().stdVector();
  const auto fine_phys_JxW = current_assembly.JxWFace().stdVector();
 
  // Find elem's own side-local reference coordinates for those same physical points.
  const Elem & elem_side_elem = _elem_side_builder(*elem, side);
  std::vector<Point> elem_side_ref_points;
  libMesh::FEMap::inverse_map(
      elem_side_elem.dim(), &elem_side_elem, fine_phys_points, elem_side_ref_points);
 
  // Measure elem's own local Jacobian at those points using dummy unit weights, then rescale so
  // that elem's face JxW, once bound for real below, exactly matches fine_neighbor's physical
  // weight computed above.
  current_assembly.reinitElemFaceRef(elem, side, TOLERANCE, &elem_side_ref_points, nullptr);
  const auto elem_dummy_JxW = current_assembly.JxWFace().stdVector();
  std::vector<Real> elem_ref_weights(fine_phys_JxW.size());
  for (const auto qp : index_range(elem_ref_weights))
    elem_ref_weights[qp] = fine_phys_JxW[qp] / elem_dummy_JxW[qp];
 
  // Bind elem's (primary side) variable to its own dofs, with the face reinitialized at the
  // corrected points/weights above. elem's residual/Jacobian contribution must accumulate
  // additively across every fine neighbor of this (elem, side) visit, since each one only covers
  // a fraction of elem's face; that local storage was already correctly prepared and zeroed
  // exactly once for this element by FEProblemBase::prepare() before this side loop ever started.
  // The loop below calls Assembly::reinitElemFaceRef() directly rather than the
  // FEProblemBase/SubProblem method of the same name: that convenience bundle also calls
  // Assembly::prepareResidual()/prepareJacobianBlock(), which would unconditionally re-zero elem's
  // local residual/Jacobian storage, discarding the prior fine neighbors' contributions as well as
  // every contribution already accumulated earlier in this element's visit (volumetric Kernels,
  // IntegratedBCs and InterfaceKernels on other sides, etc.). elem's dof indices don't change
  // across fine neighbors, so reinitializing the face at the new points/weights and recomputing
  // variable values from them is all that is needed; no further dof preparation step is required.
  for (const auto sys_num : make_range(_fe_problem.numNonlinearSystems()))
  {
    _fe_problem.assembly(_tid, sys_num)
        .reinitElemFaceRef(elem, side, TOLERANCE, &elem_side_ref_points, &elem_ref_weights);
    _fe_problem.getNonlinearSystemBase(sys_num).reinitElemFace(elem, side, _tid);
  }
  // Auxiliary variables use the system 0 assembly reinitialized above
  _fe_problem.getAuxiliarySystem().reinitElemFace(elem, side, _tid);
 
  SwapBackSentinel face_sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
  _fe_problem.reinitMaterialsFaceOnBoundary(bnd_id, elem->subdomain_id(), _tid);
  _fe_problem.reinitMaterialsBoundary(bnd_id, _tid);
 
  // Bind the neighbor variable to fine_neighbor's own dofs at the matching physical points.
  _fe_problem.setNeighborSubdomainID(fine_neighbor, _tid);
  _fe_problem.reinitNeighborPhys(fine_neighbor, fine_side, fine_phys_points, _tid);
 
  SwapBackSentinel neighbor_sentinel(_fe_problem, &FEProblem::swapBackMaterialsNeighbor, _tid);
  _fe_problem.reinitMaterialsNeighborOnBoundary(bnd_id, fine_neighbor->subdomain_id(), _tid);
 
  _fe_problem.reinitMaterialsInterface(bnd_id, _tid);
 
  computeOnInterface(bnd_id);
 
  accumulateNeighbor();
}
 
void
NonlinearThread::computeOnInterface(BoundaryID bnd_id)
{
  const auto & int_ks = _ik_warehouse->getActiveBoundaryObjects(bnd_id, _tid);
  for (const auto & interface_kernel : int_ks)
    compute(*interface_kernel);
}
 
void
NonlinearThread::onInternalSide(const Elem * elem, unsigned int side)
{
  if (_should_execute_dg)
  {
    // Pointer to the neighbor we are currently working on.
    const Elem * neighbor = elem->neighbor_ptr(side);
 
    _fe_problem.reinitElemNeighborAndLowerD(elem, side, _tid);
 
    // Set up Sentinels so that, even if one of the reinitMaterialsXXX() calls throws, we
    // still remember to swap back during stack unwinding.
    SwapBackSentinel face_sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
    _fe_problem.reinitMaterialsFace(elem->subdomain_id(), _tid);
 
    SwapBackSentinel neighbor_sentinel(_fe_problem, &FEProblem::swapBackMaterialsNeighbor, _tid);
    _fe_problem.reinitMaterialsNeighbor(neighbor->subdomain_id(), _tid);
 
    computeOnInternalFace(neighbor);
 
    accumulateNeighborLower();
  }
  if (_subdomain_has_hdg)
  {
    // Set up Sentinel class so that, after we swap in reinitMaterialsFace in prepareFace, even if
    // one of our callees throws we remember to swap back during stack unwinding. We put our
    // sentinel here as opposed to in prepareFace because we certainly don't want our materials
    // swapped back before we proceed to residual/Jacobian computation
    SwapBackSentinel sentinel(_fe_problem, &FEProblem::swapBackMaterialsFace, _tid);
 
    prepareFace(elem, side, Moose::INVALID_BOUNDARY_ID, nullptr);
    computeOnInternalFace();
  }
}
 
void
NonlinearThread::onExternalSide(const Elem * elem, unsigned int side)
{
  // Check that we don't have any interface kernels defined on this boundary
  const auto boundary_ids = _mesh.getBoundaryIDs(elem, side);
 
  bool has_interface_kernels = false;
  for (const auto bid : boundary_ids)
    if (_ik_warehouse && _ik_warehouse->hasActiveBoundaryObjects(bid, _tid))
      has_interface_kernels = true;
 
  if (has_interface_kernels)
    mooseError("Element ",
               elem->id(),
               " on side ",
               side,
               " is missing a neighbor (hence identified as an external side) but "
               "has interface kernel(s) defined on the boundary.");
}
 
void
NonlinearThread::computeOnInternalFace(const Elem * neighbor)
{
  const auto & dgks = _dg_warehouse->getActiveBlockObjects(_subdomain, _tid);
  for (const auto & dg_kernel : dgks)
    if (dg_kernel->hasBlocks(neighbor->subdomain_id()))
      compute(*dg_kernel, neighbor);
}
 
void
NonlinearThread::compute(KernelBase & kernel)
{
  compute(cast_ref<ResidualObject &>(kernel));
}
 
void
NonlinearThread::compute(FVElementalKernel & kernel)
{
  compute(cast_ref<ResidualObject &>(kernel));
}
 
void
NonlinearThread::compute(IntegratedBCBase & bc)
{
  compute(cast_ref<ResidualObject &>(bc));
}
 
void
NonlinearThread::compute(DGKernelBase & dg, const Elem * /*neighbor*/)
{
  compute(cast_ref<ResidualObject &>(dg));
}
 
void
NonlinearThread::compute(InterfaceKernelBase & ik)
{
  compute(cast_ref<ResidualObject &>(ik));
}
 
void
NonlinearThread::postElement(const Elem * /*elem*/)
{
  accumulate();
}
 
void
NonlinearThread::post()
{
  clearVarsAndMaterials();
}
 
void
NonlinearThread::printGeneralExecutionInformation() const
{
  if (_fe_problem.shouldPrintExecution(_tid))
  {
    const auto & console = _fe_problem.console();
    const auto execute_on = _fe_problem.getCurrentExecuteOnFlag();
    console << "[DBG] Beginning elemental loop to compute " + objectType() + " on " << execute_on
            << std::endl;
    mooseDoOnce(
        console << "[DBG] Execution order on each element:" << std::endl;
        console << "[DBG] - kernels on element quadrature points" << std::endl;
        console << "[DBG] - finite volume elemental kernels on element" << std::endl;
        console << "[DBG] - integrated boundary conditions on element side quadrature points"
                << std::endl;
        console << "[DBG] - DG kernels on element side quadrature points" << std::endl;
        console << "[DBG] - interface kernels on element side quadrature points" << std::endl;);
  }
}
 
void
NonlinearThread::printBlockExecutionInformation() const
{
  // Number of objects executing is approximated by size of warehouses
  const int num_objects = _kernels.size() + _fv_kernels.size() + _integrated_bcs.size() +
                          _dg_kernels.size() + _interface_kernels.size();
  const auto & console = _fe_problem.console();
  const auto block_name = _mesh.getSubdomainName(_subdomain);
 
  if (_fe_problem.shouldPrintExecution(_tid) && num_objects > 0)
  {
    if (_blocks_exec_printed.count(_subdomain))
      return;
    console << "[DBG] Ordering of " + objectType() + " Objects on block " << block_name << " ("
            << _subdomain << ")" << std::endl;
    if (_kernels.hasActiveBlockObjects(_subdomain, _tid))
    {
      console << "[DBG] Ordering of kernels:" << std::endl;
      console << _kernels.activeObjectsToFormattedString() << std::endl;
    }
    if (_fv_kernels.size())
    {
      console << "[DBG] Ordering of FV elemental kernels:" << std::endl;
      std::string fvkernels =
          std::accumulate(_fv_kernels.begin() + 1,
                          _fv_kernels.end(),
                          _fv_kernels[0]->name(),
                          [](const std::string & str_out, FVElementalKernel * kernel)
                          { return str_out + " " + kernel->name(); });
      console << ConsoleUtils::formatString(fvkernels, "[DBG]") << std::endl;
    }
    if (_dg_kernels.hasActiveBlockObjects(_subdomain, _tid))
    {
      console << "[DBG] Ordering of DG kernels:" << std::endl;
      console << _dg_kernels.activeObjectsToFormattedString() << std::endl;
    }
  }
  else if (_fe_problem.shouldPrintExecution(_tid) && num_objects == 0 &&
           !_blocks_exec_printed.count(_subdomain))
    console << "[DBG] No Active " + objectType() + " Objects on block " << block_name << " ("
            << _subdomain << ")" << std::endl;
 
  _blocks_exec_printed.insert(_subdomain);
}
 
void
NonlinearThread::printBoundaryExecutionInformation(const unsigned int bid) const
{
  if (!_fe_problem.shouldPrintExecution(_tid) || _boundaries_exec_printed.count(bid) ||
      (!_integrated_bcs.hasActiveBoundaryObjects(bid, _tid) &&
       !_interface_kernels.hasActiveBoundaryObjects(bid, _tid)))
    return;
 
  const auto & console = _fe_problem.console();
  const auto b_name = _mesh.getBoundaryName(bid);
  console << "[DBG] Ordering of " + objectType() + " Objects on boundary " << b_name << " (" << bid
          << ")" << std::endl;
 
  if (_integrated_bcs.hasActiveBoundaryObjects(bid, _tid))
  {
    console << "[DBG] Ordering of integrated boundary conditions:" << std::endl;
    console << _integrated_bcs.activeObjectsToFormattedString() << std::endl;
  }
 
  // We have not checked if we have a neighbor. This could be premature for saying we are executing
  // interface kernels. However, we should assume the execution will happen on another side of the
  // same boundary
  if (_interface_kernels.hasActiveBoundaryObjects(bid, _tid))
  {
    console << "[DBG] Ordering of interface kernels:" << std::endl;
    console << _interface_kernels.activeObjectsToFormattedString() << std::endl;
  }
 
  _boundaries_exec_printed.insert(bid);
}
 
void
NonlinearThread::prepareFace(const Elem * const elem,
                             const unsigned int side,
                             const BoundaryID bnd_id,
                             const Elem * const lower_d_elem)
{
  _fe_problem.reinitElemFace(elem, side, _tid);
 
  // Needed to use lower-dimensional variables on Materials
  if (lower_d_elem)
    _fe_problem.reinitLowerDElem(lower_d_elem, _tid);
 
  if (bnd_id != Moose::INVALID_BOUNDARY_ID)
  {
    _fe_problem.reinitMaterialsFaceOnBoundary(bnd_id, elem->subdomain_id(), _tid);
    _fe_problem.reinitMaterialsBoundary(bnd_id, _tid);
  }
  // Currently only used by HDG
  else
    _fe_problem.reinitMaterialsFace(elem->subdomain_id(), _tid);
}
 
bool
NonlinearThread::shouldComputeInternalSide(const Elem & elem, const Elem & neighbor) const
{
  // ThreadedElementLoop<ConstElemRange>::shouldComputeInternalSide gets expensive on high
  // h-refinement cases so we avoid it if possible
  _should_execute_dg = false;
  if (_subdomain_has_dg)
    _should_execute_dg =
        ThreadedElementLoop<ConstElemRange>::shouldComputeInternalSide(elem, neighbor);
  return _subdomain_has_hdg || _should_execute_dg;
}