Seperating nonlocal coupling from local one (idaholab#5913)

framework/include/base/Assembly.h

@@ -331,6 +331,9 @@ class Assembly
 
   void init();
 
+  /// Create pair of vaiables requiring nonlocal jacobian contibiutions
+  void initNonlocalCoupling();
+
   /**
    * Whether or not this assembly should utilize FE shape function caching.
    *
@@ -415,7 +418,6 @@ class Assembly
    */
   void cacheJacobianNonlocal();
 
-
   /**
    * Takes the values that are currently in the neighbor Dense Matrices and appends them to the cached values.
    */
@@ -438,6 +440,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 +561,10 @@ class Assembly
   SystemBase & _sys;
   /// Reference to coupling matrix
   CouplingMatrix * & _cm;
+  const 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

@@ -38,7 +38,6 @@
 class DisplacedProblem;
 class FEProblem;
 class MooseMesh;
-class SystemBase;
 class NonlinearSystem;
 class RandomInterface;
 class RandomData;
@@ -154,9 +153,13 @@ class FEProblem :
   void setCouplingMatrix(CouplingMatrix * cm);
   CouplingMatrix * & couplingMatrix() { return _cm; }
 
+  /// Set custom coupling matrix for variables requiring nonlocal contribution
+  void setNonlocalCouplingMatrix();
+
   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
@@ -1250,6 +1253,11 @@ class FEProblem :
   Moose::PetscSupport::PetscOptions _petsc_options;
 #endif //LIBMESH_HAVE_PETSC
 
+/**
+ * Method for sorting the MooseVariables based on variable numbers
+ */
+static bool sortMooseVariables(MooseVariable * a, MooseVariable * b) { return a->number() < b->number(); }
+
 private:
   bool _use_legacy_uo_aux_computation;
   bool _use_legacy_uo_initialization;

framework/include/base/MooseVariableBase.h

@@ -138,8 +138,10 @@ class MooseVariableBase
   const DofMap & _dof_map;
   /// DOF indices
   std::vector<dof_id_type> _dof_indices;
+
   /// mesh the variable is active in
   MooseMesh & _mesh;
+
   /// scaling factor for this variable
   Real _scaling_factor;
 };

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

@@ -21,6 +21,9 @@
 #include "GeometricSearchData.h"
 #include "MooseVariableBase.h" // VariableValue
 
+// libMesh includes
+#include "libmesh/coupling_matrix.h"
+
 class MooseMesh;
 class SubProblem;
 class Factory;
@@ -34,6 +37,7 @@ namespace libMesh
 {
 class EquationSystems;
 class DofMap;
+class CouplingMatrix;
 template <typename T> class SparseMatrix;
 template <typename T> class NumericVector;
 }
@@ -320,11 +324,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;
+  const CouplingMatrix & nonlocalCouplingMatrix() const { 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/KernelBase.h

@@ -79,6 +79,18 @@ class KernelBase :
    */
   virtual void computeOffDiagJacobianScalar(unsigned int jvar) = 0;
 
+  /**
+   * Compute this Kernel's contribution to the diagonal Jacobian entries
+   * corresponding to nonlocal dofs of the variable
+   */
+  virtual void computeNonlocalJacobian() {}
+
+  /**
+   * Computes d-residual / d-jvar... corresponding to nonlocal dofs of the jvar
+   * and stores the result in nonlocal ke
+   */
+  virtual void computeNonlocalOffDiagJacobian(unsigned int /* jvar */) {}
+
   /// Returns the variable number that this Kernel operates on.
   MooseVariable & variable();
 

framework/include/kernels/NonlocalKernel.h

@@ -22,14 +22,20 @@ class NonlocalKernel;
 template<>
 InputParameters validParams<NonlocalKernel>();
 
-class NonlocalKernel :
-  public Kernel
+/**
+ * NonlocalKernel is used for solving integral terms in integro-differential equations.
+ * Integro-differential equations includes spatial integral terms over variables in the domain.
+ * In this case the jacobian calculation is not restricted to local dofs of an element, it requires
+ * additional contributions from all the dofs in the domain. NonlocalKernel adds nonlocal jacobians
+ * to the local jacobians calculated in kernels.
+ */
+class NonlocalKernel : public Kernel
 {
 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()
+{
+  _cm_nonlocal_entry.clear();
+  const std::vector<MooseVariable *> & vars = _sys.getVariables(_tid);
+  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,26 +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();
-
-      if (vi == var->number() || vj == var->number())
-        jacobianBlockNonlocal(vi,vj).resize(ivar.dofIndices().size(), jvar.allDofIndices().size());
-    }
-  }
+    unsigned int vi = ivar.number();
+    unsigned int vj = jvar.number();
 
-  for (unsigned int i = 0; i < _sub_Re.size(); i++)
-  {
-    _sub_Re[i][var->number()].resize(var->dofIndices().size());
-    _sub_Re[i][var->number()].zero();
+    if (vi == var->number() || vj == var->number())
+      jacobianBlockNonlocal(vi,vj).resize(ivar.dofIndices().size(), jvar.allDofIndices().size());
   }
 }
 
@@ -1452,8 +1456,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 +1506,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 +1582,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();
       }
   }
 }