DisplacedProblem.h

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

#pragma once

#include "SubProblem.h"
#include "DisplacedSystem.h"
#include "GeometricSearchData.h"
#include "ThreadedNodeLoop.h"

// libMesh
#include "libmesh/equation_systems.h"
#include "libmesh/enum_quadrature_type.h"

// Forward declarations
class MooseVariableFieldBase;
class AssemblyData;
class MooseMesh;
class Assembly;
class FEProblemBase;
class LineSearch;

// libMesh forward declarations
namespace libMesh
{
template <typename T>
class NumericVector;
class DofMap;
}

class DisplacedProblem : public SubProblem
{
public:
  static InputParameters validParams();

  DisplacedProblem(DisplacedProblem &&) = delete;
  DisplacedProblem & operator=(DisplacedProblem &&) = delete;

  DisplacedProblem(const InputParameters & parameters);
  ~DisplacedProblem();

  virtual EquationSystems & es() override { return _eq; }
  virtual MooseMesh & mesh() override { return _mesh; }
  virtual const MooseMesh & mesh() const override { return _mesh; }
  const MooseMesh & mesh(bool libmesh_dbg_var(use_displaced)) const override
  {
    mooseAssert(use_displaced, "An undisplaced mesh was queried from the displaced problem");
    return mesh();
  }
  MooseMesh & refMesh();

  DisplacedSystem & solverSys(const unsigned int sys_num);
  DisplacedSystem & auxSys() { return *_displaced_aux; }

  virtual const SystemBase & systemBaseNonlinear(const unsigned int sys_num) const override;
  virtual SystemBase & systemBaseNonlinear(const unsigned int sys_num) override;

  virtual const SystemBase & systemBaseLinear(const unsigned int sys_num) const override;
  virtual SystemBase & systemBaseLinear(const unsigned int sys_num) override;

  virtual const SystemBase & systemBaseSolver(const unsigned int sys_num) const override;
  virtual SystemBase & systemBaseSolver(const unsigned int sys_num) override;

  virtual const SystemBase & systemBaseAuxiliary() const override { return *_displaced_aux; }
  virtual SystemBase & systemBaseAuxiliary() override { return *_displaced_aux; }

  // Return a constant reference to the vector of variable names.
  const std::vector<std::string> & getDisplacementVarNames() const { return _displacements; }

  virtual void createQRules(QuadratureType type,
                            Order order,
                            Order volume_order,
                            Order face_order,
                            SubdomainID block,
                            bool allow_negative_qweights = true);

  void bumpVolumeQRuleOrder(Order order, SubdomainID block);
  void bumpAllQRuleOrder(Order order, SubdomainID block);

  virtual void init() override;
  virtual bool solverSystemConverged(const unsigned int solver_sys_num) override;
  virtual unsigned int nlSysNum(const NonlinearSystemName & nl_sys_name) const override;
  virtual unsigned int linearSysNum(const LinearSystemName & sys_name) const override;
  virtual unsigned int solverSysNum(const SolverSystemName & sys_name) const override;

  void addTimeIntegrator();

  virtual void saveOldSolutions();

  virtual void restoreOldSolutions();

  void syncAuxSolution(const NumericVector<Number> & aux_soln);

  void syncSolutions();

  void syncSolutions(const std::map<unsigned int, const NumericVector<Number> *> & nl_solns,
                     const NumericVector<Number> & aux_soln);

  virtual void updateMesh(bool mesh_changing = false);

  virtual void updateMesh(const std::map<unsigned int, const NumericVector<Number> *> & nl_soln,
                          const NumericVector<Number> & aux_soln);

  virtual TagID addVectorTag(const TagName & tag_name,
                             const Moose::VectorTagType type = Moose::VECTOR_TAG_RESIDUAL) override;
  virtual const VectorTag & getVectorTag(const TagID tag_id) const override;
  virtual TagID getVectorTagID(const TagName & tag_name) const override;
  virtual TagName vectorTagName(const TagID tag_id) const override;
  virtual bool vectorTagExists(const TagID tag_id) const override;
  virtual bool vectorTagExists(const TagName & tag_name) const override;
  virtual unsigned int
  numVectorTags(const Moose::VectorTagType type = Moose::VECTOR_TAG_ANY) const override;
  virtual const std::vector<VectorTag> &
  getVectorTags(const Moose::VectorTagType type = Moose::VECTOR_TAG_ANY) const override;
  virtual Moose::VectorTagType vectorTagType(const TagID tag_id) const override;

  virtual TagID addMatrixTag(TagName tag_name) override;
  virtual TagID getMatrixTagID(const TagName & tag_name) const override;
  virtual TagName matrixTagName(TagID tag) override;
  virtual bool matrixTagExists(const TagName & tag_name) const override;
  virtual bool matrixTagExists(TagID tag_id) const override;
  virtual unsigned int numMatrixTags() const override;
  virtual bool safeAccessTaggedMatrices() const override;
  virtual bool safeAccessTaggedVectors() const override;

  virtual bool isTransient() const override;

  // Variables /////
  virtual bool hasVariable(const std::string & var_name) const override;
  using SubProblem::getVariable;
  virtual const MooseVariableFieldBase &
  getVariable(const THREAD_ID tid,
              const std::string & var_name,
              Moose::VarKindType expected_var_type = Moose::VarKindType::VAR_ANY,
              Moose::VarFieldType expected_var_field_type =
                  Moose::VarFieldType::VAR_FIELD_ANY) const override;
  virtual MooseVariable & getStandardVariable(const THREAD_ID tid,
                                              const std::string & var_name) override;
  virtual MooseVariableFieldBase & getActualFieldVariable(const THREAD_ID tid,
                                                          const std::string & var_name) override;
  virtual VectorMooseVariable & getVectorVariable(const THREAD_ID tid,
                                                  const std::string & var_name) override;
  virtual ArrayMooseVariable & getArrayVariable(const THREAD_ID tid,
                                                const std::string & var_name) override;
  virtual bool hasScalarVariable(const std::string & var_name) const override;
  virtual MooseVariableScalar & getScalarVariable(const THREAD_ID tid,
                                                  const std::string & var_name) override;
  virtual System & getSystem(const std::string & var_name) override;

  virtual void addVariable(const std::string & var_type,
                           const std::string & name,
                           InputParameters & parameters,
                           unsigned int nl_system_number);
  virtual void addAuxVariable(const std::string & var_type,
                              const std::string & name,
                              InputParameters & parameters);
  //
  // Adaptivity /////
  virtual void initAdaptivity();
  void meshChanged(bool contract_mesh, bool clean_refinement_flags);

  // reinit /////
  virtual void prepare(const Elem * elem, const THREAD_ID tid) override;
  virtual void prepareNonlocal(const THREAD_ID tid);
  virtual void prepareFace(const Elem * elem, const THREAD_ID tid) override;
  virtual void prepare(const Elem * elem,
                       unsigned int ivar,
                       unsigned int jvar,
                       const std::vector<dof_id_type> & dof_indices,
                       const THREAD_ID tid) override;
  virtual void setCurrentSubdomainID(const Elem * elem, const THREAD_ID tid) override;
  virtual void
  setNeighborSubdomainID(const Elem * elem, unsigned int side, const THREAD_ID tid) override;
  virtual void prepareBlockNonlocal(unsigned int ivar,
                                    unsigned int jvar,
                                    const std::vector<dof_id_type> & idof_indices,
                                    const std::vector<dof_id_type> & jdof_indices,
                                    const THREAD_ID tid);
  virtual void prepareAssembly(const THREAD_ID tid) override;
  virtual void prepareAssemblyNeighbor(const THREAD_ID tid);

  virtual bool reinitDirac(const Elem * elem, const THREAD_ID tid) override;

  virtual void reinitElem(const Elem * elem, const THREAD_ID tid) override;
  virtual void reinitElemPhys(const Elem * elem,
                              const std::vector<Point> & phys_points_in_elem,
                              const THREAD_ID tid) override;
  virtual void reinitElemFace(const Elem * elem, unsigned int side, const THREAD_ID tid) override;
  virtual void reinitNode(const Node * node, const THREAD_ID tid) override;
  virtual void reinitNodeFace(const Node * node, BoundaryID bnd_id, const THREAD_ID tid) override;
  virtual void reinitNodes(const std::vector<dof_id_type> & nodes, const THREAD_ID tid) override;
  virtual void reinitNodesNeighbor(const std::vector<dof_id_type> & nodes,
                                   const THREAD_ID tid) override;
  virtual void reinitNeighbor(const Elem * elem, unsigned int side, const THREAD_ID tid) override;

  void reinitNeighbor(const Elem * elem,
                      unsigned int side,
                      const THREAD_ID tid,
                      const std::vector<Point> * neighbor_reference_points);

  virtual void reinitNeighborPhys(const Elem * neighbor,
                                  unsigned int neighbor_side,
                                  const std::vector<Point> & physical_points,
                                  const THREAD_ID tid) override;
  virtual void reinitNeighborPhys(const Elem * neighbor,
                                  const std::vector<Point> & physical_points,
                                  const THREAD_ID tid) override;
  virtual void
  reinitElemNeighborAndLowerD(const Elem * elem, unsigned int side, const THREAD_ID tid) override;
  virtual void reinitScalars(const THREAD_ID tid,
                             bool reinit_for_derivative_reordering = false) override;
  virtual void reinitOffDiagScalars(const THREAD_ID tid) override;

  virtual void getDiracElements(std::set<const Elem *> & elems) override;
  virtual void clearDiracInfo() override;

  virtual void addResidual(const THREAD_ID tid) override;
  virtual void addResidualNeighbor(const THREAD_ID tid) override;
  virtual void addResidualLower(const THREAD_ID tid) override;

  virtual void addCachedResidualDirectly(NumericVector<Number> & residual, const THREAD_ID tid);

  virtual void setResidual(NumericVector<Number> & residual, const THREAD_ID tid) override;
  virtual void setResidualNeighbor(NumericVector<Number> & residual, const THREAD_ID tid) override;

  virtual void addJacobian(const THREAD_ID tid) override;
  virtual void addJacobianNonlocal(const THREAD_ID tid);
  virtual void addJacobianNeighbor(const THREAD_ID tid) override;
  virtual void addJacobianNeighborLowerD(const THREAD_ID tid) override;
  virtual void addJacobianLowerD(const THREAD_ID tid) override;
  virtual void addJacobianBlockTags(SparseMatrix<Number> & jacobian,
                                    unsigned int ivar,
                                    unsigned int jvar,
                                    const libMesh::DofMap & dof_map,
                                    std::vector<dof_id_type> & dof_indices,
                                    const std::set<TagID> & tags,
                                    const THREAD_ID tid);
  void addJacobianBlockNonlocal(SparseMatrix<Number> & jacobian,
                                unsigned int ivar,
                                unsigned int jvar,
                                const libMesh::DofMap & dof_map,
                                const std::vector<dof_id_type> & idof_indices,
                                const std::vector<dof_id_type> & jdof_indices,
                                const std::set<TagID> & tags,
                                const THREAD_ID tid);
  virtual void addJacobianNeighbor(SparseMatrix<Number> & jacobian,
                                   unsigned int ivar,
                                   unsigned int jvar,
                                   const libMesh::DofMap & dof_map,
                                   std::vector<dof_id_type> & dof_indices,
                                   std::vector<dof_id_type> & neighbor_dof_indices,
                                   const std::set<TagID> & tags,
                                   const THREAD_ID tid) override;

  virtual void cacheJacobianNonlocal(const THREAD_ID tid);

  virtual void prepareShapes(unsigned int var, const THREAD_ID tid) override;
  virtual void prepareFaceShapes(unsigned int var, const THREAD_ID tid) override;
  virtual void prepareNeighborShapes(unsigned int var, const THREAD_ID tid) override;

  virtual Assembly & assembly(const THREAD_ID tid, const unsigned int sys_num) override;
  virtual const Assembly & assembly(const THREAD_ID tid, const unsigned int sys_num) const override;

  // Geom Search /////
  virtual void updateGeomSearch(
      GeometricSearchData::GeometricSearchType type = GeometricSearchData::ALL) override;
  virtual GeometricSearchData & geomSearchData() override { return _geometric_search_data; }

  virtual bool computingPreSMOResidual(const unsigned int nl_sys_num) const override;

  virtual void onTimestepBegin() override;
  virtual void onTimestepEnd() override;

  virtual std::set<dof_id_type> & ghostedElems() override;

  virtual void addGhostedElem(dof_id_type elem_id) override;

  virtual void addGhostedBoundary(BoundaryID boundary_id) override;

  virtual void ghostGhostedBoundaries() override;

  void undisplaceMesh();

  virtual LineSearch * getLineSearch() override;

  virtual const CouplingMatrix * couplingMatrix(const unsigned int nl_sys_num) const override;

  virtual bool haveDisplaced() const override final { return true; }

  virtual bool computingScalingJacobian() const override final;

  virtual bool computingScalingResidual() const override final;

  virtual void initialSetup() override;
  virtual void timestepSetup() override;
  virtual void customSetup(const ExecFlagType & exec_type) override;
  virtual void residualSetup() override;
  virtual void jacobianSetup() override;

  using SubProblem::haveADObjects;
  virtual void haveADObjects(bool have_ad_objects) override;

  virtual std::size_t numNonlinearSystems() const override;

  virtual std::size_t numLinearSystems() const override;

  virtual std::size_t numSolverSystems() const override;

  virtual unsigned int currentNlSysNum() const override;

  virtual unsigned int currentLinearSysNum() const override;

  virtual const std::vector<VectorTag> & currentResidualVectorTags() const override;

  virtual void needFV() override;
  virtual bool haveFV() const override;

  virtual bool hasNonlocalCoupling() const override;
  virtual bool checkNonlocalCouplingRequirement() const override;
  virtual const libMesh::CouplingMatrix & nonlocalCouplingMatrix(const unsigned i) const override;

protected:
  FEProblemBase & _mproblem;
  MooseMesh & _mesh;
  EquationSystems _eq;
  MooseMesh & _ref_mesh;
  std::vector<std::string> _displacements;

  std::vector<std::unique_ptr<DisplacedSystem>> _displaced_solver_systems;
  std::unique_ptr<DisplacedSystem> _displaced_aux;

  std::vector<const NumericVector<Number> *> _nl_solution;

  const NumericVector<Number> * _aux_solution;

  std::vector<std::vector<std::unique_ptr<Assembly>>> _assembly;

  GeometricSearchData _geometric_search_data;

  class UpdateDisplacedMeshThread : public ThreadedNodeLoop<NodeRange, NodeRange::const_iterator>
  {
  public:
    UpdateDisplacedMeshThread(FEProblemBase & fe_problem, DisplacedProblem & displaced_problem);

    UpdateDisplacedMeshThread(UpdateDisplacedMeshThread & x, Threads::split split);

    virtual void onNode(NodeRange::const_iterator & nd) override;

    void join(const UpdateDisplacedMeshThread & y)
    {
      if (y._has_displacement)
        _has_displacement = true;
    }

    bool hasDisplacement()
    {
      mooseAssert(!Threads::in_threads,
                  "This function requires a MPI all-gathering operation that cannot be in a "
                  "threaded scope.");
      _ref_mesh.comm().max(_has_displacement);
      return _has_displacement;
    }

  protected:
    void init();

    DisplacedProblem & _displaced_problem;
    MooseMesh & _ref_mesh;
    const std::vector<const NumericVector<Number> *> & _nl_soln;
    const NumericVector<Number> & _aux_soln;

    // Solution vectors with expanded ghosting, for ReplicatedMesh or
    // for DistributedMesh cases where the standard algebraic ghosting
    // doesn't reach as far as the geometric ghosting
    std::map<unsigned int,
             std::pair<const NumericVector<Number> *, std::shared_ptr<NumericVector<Number>>>>
        _sys_to_nonghost_and_ghost_soln;

  private:
    std::map<unsigned int, std::pair<std::vector<unsigned int>, std::vector<unsigned int>>>
        _sys_to_var_num_and_direction;

    bool _has_displacement;
  };

private:
  virtual std::pair<bool, unsigned int>
  determineSolverSystem(const std::string & var_name,
                        bool error_if_not_found = false) const override;

  friend class UpdateDisplacedMeshThread;
  friend class Restartable;
};

inline DisplacedSystem &
DisplacedProblem::solverSys(const unsigned int sys_num)
{
  mooseAssert(sys_num < _displaced_solver_systems.size(),
              "System number greater than the number of nonlinear systems");
  return *_displaced_solver_systems[sys_num];
}

inline const SystemBase &
DisplacedProblem::systemBaseNonlinear(const unsigned int sys_num) const
{
  mooseAssert(sys_num < _displaced_solver_systems.size(),
              "System number greater than the number of nonlinear systems");
  return *_displaced_solver_systems[sys_num];
}

inline SystemBase &
DisplacedProblem::systemBaseNonlinear(const unsigned int sys_num)
{
  mooseAssert(sys_num < _displaced_solver_systems.size(),
              "System number greater than the number of nonlinear systems");
  return *_displaced_solver_systems[sys_num];
}

inline const SystemBase &
DisplacedProblem::systemBaseLinear(const unsigned int /*sys_num*/) const
{
  mooseError("Linear systems are not supported for displaced problems yet.");
}

inline SystemBase &
DisplacedProblem::systemBaseLinear(const unsigned int /*sys_num*/)
{
  mooseError("Linear systems are not supported for displaced problems yet.");
}

inline const SystemBase &
DisplacedProblem::systemBaseSolver(const unsigned int sys_num) const
{
  mooseAssert(sys_num < _displaced_solver_systems.size(),
              "System number greater than the number of solver systems");
  return *_displaced_solver_systems[sys_num];
}

inline SystemBase &
DisplacedProblem::systemBaseSolver(const unsigned int sys_num)
{
  mooseAssert(sys_num < _displaced_solver_systems.size(),
              "System number greater than the number of solver systems");
  return *_displaced_solver_systems[sys_num];
}