Seperating nonlocal coupling from local one (idaholab#5913)

framework/include/base/Assembly.h

@@ -330,6 +330,7 @@ class Assembly
   void reinitNodeNeighbor(const Node * node);
 
   void init();
+  void initNonlocalCoupling();
 
   /**
    * Whether or not this assembly should utilize FE shape function caching.
@@ -438,6 +439,7 @@ class Assembly
   void cacheJacobianBlockNonlocal(DenseMatrix<Number> & jac_block, const std::vector<dof_id_type> & idof_indices, const std::vector<dof_id_type> & jdof_indices, Real scaling_factor);
 
   std::vector<std::pair<MooseVariable *, MooseVariable *> > & couplingEntries() { return _cm_entry; }
+  std::vector<std::pair<MooseVariable *, MooseVariable *> > & nonlocalCouplingEntries() { return _cm_nonlocal_entry; }
 
   const VariablePhiValue & phi() { return _phi; }
   const VariablePhiGradient & gradPhi() { return _grad_phi; }
@@ -558,8 +560,10 @@ class Assembly
   SystemBase & _sys;
   /// Reference to coupling matrix
   CouplingMatrix * & _cm;
+  CouplingMatrix * & _nonlocal_cm;
   /// Entries in the coupling matrix (only for field variables)
   std::vector<std::pair<MooseVariable *, MooseVariable *> > _cm_entry;
+  std::vector<std::pair<MooseVariable *, MooseVariable *> > _cm_nonlocal_entry;
   /// Flag that indicates if the jacobian block was used
   std::vector<std::vector<unsigned char> > _jacobian_block_used;
   std::vector<std::vector<unsigned char> > _jacobian_block_nonlocal_used;

framework/include/base/FEProblem.h

@@ -154,9 +154,13 @@ class FEProblem :
   void setCouplingMatrix(CouplingMatrix * cm);
   CouplingMatrix * & couplingMatrix() { return _cm; }
 
+  void setNonlocalCouplingMatrix(CouplingMatrix * cm);
+  CouplingMatrix * & nonlocalCouplingMatrix() { return _nonlocal_cm; }
+
   bool areCoupled(unsigned int ivar, unsigned int jvar);
 
   std::vector<std::pair<MooseVariable *, MooseVariable *> > & couplingEntries(THREAD_ID tid);
+  std::vector<std::pair<MooseVariable *, MooseVariable *> > & nonlocalCouplingEntries(THREAD_ID tid);
 
   /**
    * Check for converence of the nonlinear solution
@@ -1065,7 +1069,6 @@ class FEProblem :
 
   Moose::CouplingType _coupling;                        ///< Type of variable coupling
   CouplingMatrix * _cm;                                 ///< Coupling matrix for variables. It is diagonal, since we do only block diagonal preconditioning.
-
   // Dimension of the subspace spanned by the vectors with a given prefix
   std::map<std::string,unsigned int> _subspace_dim;
 

framework/include/base/NonlinearSystem.h

@@ -230,7 +230,9 @@ class NonlinearSystem : public SystemTempl<TransientNonlinearImplicitSystem>,
    */
   void constraintJacobians(SparseMatrix<Number> & jacobian, bool displaced);
 
-  const std::vector<dof_id_type> & getVariableGlobalDoFs(const std::string & var_name);
+  /// set all the global dof indices for a nonlinear variable
+  void setVariableGlobalDoFs(const std::string & var_name);
+  const std::vector<dof_id_type> & getVariableGlobalDoFs() { return _var_all_dof_indices; }
 
   /**
    * Computes Jacobian

framework/include/base/SubProblem.h

@@ -34,6 +34,7 @@ namespace libMesh
 {
 class EquationSystems;
 class DofMap;
+class CouplingMatrix;
 template <typename T> class SparseMatrix;
 template <typename T> class NumericVector;
 }
@@ -320,11 +321,14 @@ class SubProblem : public Problem
   virtual void registerRestartableData(std::string name, RestartableDataValue * data, THREAD_ID tid);
 
   std::map<std::string, std::vector<dof_id_type> > _var_dof_map;
+  CouplingMatrix * & nonlocalCouplingMatrix() { return _nonlocal_cm; }
 
 protected:
   /// The Factory for building objects
   Factory & _factory;
 
+  CouplingMatrix * _nonlocal_cm; /// nonlocal coupling matrix;
+
   /// Type of coordinate system per subdomain
   std::map<SubdomainID, Moose::CoordinateSystemType> _coord_sys;
 

framework/include/kernels/EigenKernel.h

@@ -37,6 +37,8 @@ class EigenKernel : public KernelBase
   virtual void computeJacobian() override;
   virtual void computeOffDiagJacobian(unsigned int /*jvar*/) override {}
   virtual void computeOffDiagJacobianScalar(unsigned int /*jvar*/) override {}
+  virtual void computeNonlocalJacobian() override {}
+  virtual void computeNonlocalOffDiagJacobian(unsigned int /*jvar*/) override {}
 
   EigenKernel(const InputParameters & parameters);
   virtual bool enabled() override;

framework/include/kernels/Kernel.h

@@ -33,6 +33,9 @@ class Kernel :
   virtual void computeOffDiagJacobian(unsigned int jvar) override;
   virtual void computeOffDiagJacobianScalar(unsigned int jvar) override;
 
+  virtual void computeNonlocalJacobian() override {}
+  virtual void computeNonlocalOffDiagJacobian(unsigned int /* jvar */) override {}
+
 protected:
   /// Compute this Kernel's contribution to the residual at the current quadrature point
   virtual Real computeQpResidual() = 0;

framework/include/kernels/KernelBase.h

@@ -73,6 +73,12 @@ class KernelBase :
   /// Computes d-residual / d-jvar... storing the result in Ke.
   virtual void computeOffDiagJacobian(unsigned int jvar) = 0;
 
+  /// Compute this Kernel's contribution to the diagonal Jacobian entries
+  virtual void computeNonlocalJacobian() = 0;
+
+  /// Computes d-residual / d-jvar... storing the result in Ke.
+  virtual void computeNonlocalOffDiagJacobian(unsigned int jvar) = 0;
+
   /**
    * Computes jacobian block with respect to a scalar variable
    * @param jvar The number of the scalar variable

framework/include/kernels/NonlocalKernel.h

@@ -28,8 +28,8 @@ class NonlocalKernel :
 public:
   NonlocalKernel(const InputParameters & parameters);
 
-  virtual void computeJacobian();
-  virtual void computeOffDiagJacobian(unsigned int jvar);
+  virtual void computeNonlocalJacobian();
+  virtual void computeNonlocalOffDiagJacobian(unsigned int jvar);
 
 protected:
   /// Compute this Kernel's contribution to the Jacobian corresponding to nolocal dof at the current quadrature point

framework/src/base/Assembly.C

@@ -37,6 +37,7 @@
 Assembly::Assembly(SystemBase & sys, CouplingMatrix * & cm, THREAD_ID tid) :
     _sys(sys),
     _cm(cm),
+    _nonlocal_cm(_sys.subproblem().nonlocalCouplingMatrix()),
     _dof_map(_sys.dofMap()),
     _tid(tid),
     _mesh(sys.mesh()),
@@ -955,6 +956,23 @@ Assembly::init()
   }
 }
 
+void
+Assembly::initNonlocalCoupling()
+{
+  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
+  _cm_nonlocal_entry.clear();
+  for (const auto & jvar : vars)
+  {
+    unsigned int j = jvar->number();
+    for (const auto & ivar : vars)
+    {
+      unsigned int i = ivar->number();
+      if ((*_nonlocal_cm)(i, j) != 0)
+        _cm_nonlocal_entry.push_back(std::make_pair(ivar, jvar));
+    }
+  }
+}
+
 void
 Assembly::prepare()
 {
@@ -983,21 +1001,17 @@ Assembly::prepare()
 void
 Assembly::prepareNonlocal()
 {
-  for (const auto & it : _cm_entry)
+  for (const auto & it : _cm_nonlocal_entry)
   {
     MooseVariable & ivar = *(it.first);
     MooseVariable & jvar = *(it.second);
 
-    const auto map_it = _sys.subproblem()._var_dof_map.find(jvar.name());
-    if (map_it != _sys.subproblem()._var_dof_map.end())
-    {
-      unsigned int vi = ivar.number();
-      unsigned int vj = jvar.number();
+    unsigned int vi = ivar.number();
+    unsigned int vj = jvar.number();
 
-      jacobianBlockNonlocal(vi,vj).resize(ivar.dofIndices().size(), jvar.allDofIndices().size());
-      jacobianBlockNonlocal(vi,vj).zero();
-      _jacobian_block_nonlocal_used[vi][vj] = 0;
-    }
+    jacobianBlockNonlocal(vi,vj).resize(ivar.dofIndices().size(), jvar.allDofIndices().size());
+    jacobianBlockNonlocal(vi,vj).zero();
+    _jacobian_block_nonlocal_used[vi][vj] = 0;
   }
 }
 
@@ -1026,20 +1040,16 @@ Assembly::prepareVariable(MooseVariable * var)
 void
 Assembly::prepareVariableNonlocal(MooseVariable * var)
 {
-  for (const auto & it : _cm_entry)
+  for (const auto & it : _cm_nonlocal_entry)
   {
     MooseVariable & ivar = *(it.first);
     MooseVariable & jvar = *(it.second);
 
-    const auto map_it = _sys.subproblem()._var_dof_map.find(jvar.name());
-    if (map_it != _sys.subproblem()._var_dof_map.end())
-    {
-      unsigned int vi = ivar.number();
-      unsigned int vj = jvar.number();
+    unsigned int vi = ivar.number();
+    unsigned int vj = jvar.number();
 
-      if (vi == var->number() || vj == var->number())
-        jacobianBlockNonlocal(vi,vj).resize(ivar.dofIndices().size(), jvar.allDofIndices().size());
-    }
+    if (vi == var->number() || vj == var->number())
+      jacobianBlockNonlocal(vi,vj).resize(ivar.dofIndices().size(), jvar.allDofIndices().size());
   }
 
   for (unsigned int i = 0; i < _sub_Re.size(); i++)
@@ -1452,8 +1462,9 @@ Assembly::cacheJacobianBlockNonlocal(DenseMatrix<Number> & jac_block, const std:
 void
 Assembly::addCachedJacobian(SparseMatrix<Number> & jacobian)
 {
-  mooseAssert(_cached_jacobian_rows.size() == _cached_jacobian_cols.size(),
-              "Error: Cached data sizes MUST be the same!");
+  if (!_sys.subproblem().checkNonlocalCouplingRequirement())
+    mooseAssert(_cached_jacobian_rows.size() == _cached_jacobian_cols.size(),
+                "Error: Cached data sizes MUST be the same!");
 
   for (unsigned int i=0; i<_cached_jacobian_rows.size(); i++)
     jacobian.add(_cached_jacobian_rows[i], _cached_jacobian_cols[i], _cached_jacobian_values[i]);
@@ -1501,18 +1512,15 @@ Assembly::addJacobian(SparseMatrix<Number> & jacobian)
     }
   }
 }
+
 void
 Assembly::addJacobianNonlocal(SparseMatrix<Number> & jacobian)
 {
   const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
   for (const auto & ivar : vars)
     for (const auto & jvar : vars)
-    {
-      const auto it = _sys.subproblem()._var_dof_map.find(jvar->name());
-      if (it != _sys.subproblem()._var_dof_map.end())
-        if ((*_cm)(ivar->number(), jvar->number()) != 0 && _jacobian_block_nonlocal_used[ivar->number()][jvar->number()])
-          addJacobianBlock(jacobian, jacobianBlockNonlocal(ivar->number(), jvar->number()), ivar->dofIndices(), jvar->allDofIndices(), ivar->scalingFactor());
-    }
+      if ((*_nonlocal_cm)(ivar->number(), jvar->number()) != 0 && _jacobian_block_nonlocal_used[ivar->number()][jvar->number()])
+        addJacobianBlock(jacobian, jacobianBlockNonlocal(ivar->number(), jvar->number()), ivar->dofIndices(), jvar->allDofIndices(), ivar->scalingFactor());
 }
 
 void
@@ -1580,12 +1588,8 @@ Assembly::cacheJacobianNonlocal()
   const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
   for (const auto & ivar : vars)
     for (const auto & jvar : vars)
-    {
-      const auto it = _sys.subproblem()._var_dof_map.find(jvar->name());
-      if (it != _sys.subproblem()._var_dof_map.end())
-        if ((*_cm)(ivar->number(), jvar->number()) != 0 && _jacobian_block_nonlocal_used[ivar->number()][jvar->number()])
-          cacheJacobianBlockNonlocal(jacobianBlockNonlocal(ivar->number(), jvar->number()), ivar->dofIndices(), jvar->allDofIndices(), ivar->scalingFactor());
-    }
+      if ((*_nonlocal_cm)(ivar->number(), jvar->number()) != 0 && _jacobian_block_nonlocal_used[ivar->number()][jvar->number()])
+        cacheJacobianBlockNonlocal(jacobianBlockNonlocal(ivar->number(), jvar->number()), ivar->dofIndices(), jvar->allDofIndices(), ivar->scalingFactor());
 }
 
 void

framework/src/base/ComputeFullJacobianThread.C

@@ -19,6 +19,7 @@
 #include "IntegratedBC.h"
 #include "DGKernel.h"
 #include "InterfaceKernel.h"
+#include "NonlocalKernel.h"
 // libmesh includes
 #include "libmesh/threads.h"
 
@@ -62,12 +63,39 @@ ComputeFullJacobianThread::computeJacobian()
       // only if there are dofs for j-variable (if it is subdomain restricted var, there may not be any)
       const std::vector<MooseSharedPointer<KernelBase> > & kernels = _kernels.getActiveVariableBlockObjects(ivar, _subdomain, _tid);
       for (const auto & kernel : kernels)
-      {
         if ((kernel->variable().number() == ivar) && kernel->isImplicit())
         {
           kernel->subProblem().prepareShapes(jvar, _tid);
           kernel->computeOffDiagJacobian(jvar);
         }
+    }
+  }
+
+  /// done only when nonlocal kernels exist in the system
+  if (_fe_problem.checkNonlocalCouplingRequirement())
+  {
+    std::vector<std::pair<MooseVariable *, MooseVariable *> > & cne = _fe_problem.nonlocalCouplingEntries(_tid);
+    for (const auto & it : cne)
+    {
+      MooseVariable & ivariable = *(it.first);
+      MooseVariable & jvariable = *(it.second);
+
+      unsigned int ivar = ivariable.number();
+      unsigned int jvar = jvariable.number();
+
+      if (ivariable.activeOnSubdomain(_subdomain) && jvariable.activeOnSubdomain(_subdomain) && _kernels.hasActiveVariableBlockObjects(ivar, _subdomain, _tid))
+      {
+        const std::vector<MooseSharedPointer<KernelBase> > & kernels = _kernels.getActiveVariableBlockObjects(ivar, _subdomain, _tid);
+        for (const auto & kernel : kernels)
+        {
+          MooseSharedPointer<NonlocalKernel> nonlocal_kernel = MooseSharedNamespace::dynamic_pointer_cast<NonlocalKernel>(kernel);
+          if (nonlocal_kernel)
+            if ((kernel->variable().number() == ivar) && kernel->isImplicit())
+            {
+              kernel->subProblem().prepareShapes(jvar, _tid);
+              kernel->computeNonlocalOffDiagJacobian(jvar);
+            }
+        }
       }
     }
   }

framework/src/base/ComputeJacobianThread.C

@@ -20,6 +20,7 @@
 #include "DGKernel.h"
 #include "InterfaceKernel.h"
 #include "KernelWarehouse.h"
+#include "NonlocalKernel.h"
 
 // libmesh includes
 #include "libmesh/threads.h"
@@ -64,6 +65,9 @@ ComputeJacobianThread::computeJacobian()
       {
         kernel->subProblem().prepareShapes(kernel->variable().number(), _tid);
         kernel->computeJacobian();
+        MooseSharedPointer<NonlocalKernel> nonlocal_kernel = MooseSharedNamespace::dynamic_pointer_cast<NonlocalKernel>(kernel);
+        if (nonlocal_kernel)
+          kernel->computeNonlocalJacobian();
       }
   }
 }