Merge pull request #1 from lindsayad/mortar-scalar-assembly

Make residual-and-Jacobian separate work and test

framework/src/kernels/ADKernelScalarBase.C

@@ -41,7 +41,9 @@ ADKernelScalarBase::ADKernelScalarBase(const InputParameters & parameters)
     _k_order(_use_scalar ? _kappa_var_ptr->order() : 0),
     _kappa(_use_scalar ? _kappa_var_ptr->adSln() : _kappa_dummy)
 {
-  // add some error checks here
+#ifndef MOOSE_GLOBAL_AD_INDEXING
+  mooseError("ADKernelScalarBase only supported for global AD indexing");
+#endif
 }
 
 void
@@ -74,44 +76,21 @@ ADKernelScalarBase::computeJacobian()
 
   if (_use_scalar)
   {
-#ifndef MOOSE_GLOBAL_AD_INDEXING
-    mooseError("Jacobian assembly not coded for non-default AD");
-#endif
-#ifndef MOOSE_SPARSE_AD
-    mooseError("Jacobian assembly not coded for non-sparse AD");
-#else
     computeScalarResidualsForJacobian();
     _assembly.processResidualsAndJacobian(_scalar_residuals,
                                           _kappa_var_ptr->dofIndices(),
                                           _vector_tags,
                                           _matrix_tags,
                                           _kappa_var_ptr->scalingFactor());
-    // const std::vector<std::pair<MooseVariableScalar *, MooseVariableScalar *>> kvar_kvar_coupling
-    // =
-    //     {std::make_pair(_kappa_var_ptr, _kappa_var_ptr)};
-    // computeADScalarJacobian(kvar_kvar_coupling);
-#endif
   }
 }
 
 void
 ADKernelScalarBase::computeOffDiagJacobian(const unsigned int jvar_num)
 {
-  if (_use_scalar)
-  {
-#ifndef MOOSE_GLOBAL_AD_INDEXING
-    mooseError("off-diagonal Jacobian assembly not coded for non-default AD");
-#endif
-    ADKernel::computeResidualsForJacobian();
-#ifdef MOOSE_SPARSE_AD
-    _assembly.processResidualsAndJacobian(
-        _residuals, _var.dofIndices(), _vector_tags, _matrix_tags, _var.scalingFactor());
-#endif
-  }
-  else
-  {
-    ADKernel::computeOffDiagJacobian(jvar_num); // d-_var-residual / d-jvar
-  }
+  // Only need to do this once because AD does all the derivatives at once
+  if (jvar_num == _var.number())
+    computeJacobian();
 }
 
 void
@@ -122,23 +101,17 @@ ADKernelScalarBase::computeOffDiagJacobianScalar(const unsigned int /*jvar_num*/
 void
 ADKernelScalarBase::computeResidualAndJacobian()
 {
-#ifdef MOOSE_GLOBAL_AD_INDEXING
   ADKernel::computeResidualAndJacobian();
 
   if (_use_scalar)
   {
-#ifdef MOOSE_SPARSE_AD
     computeScalarResidualsForJacobian();
     _assembly.processResidualsAndJacobian(_scalar_residuals,
                                           _kappa_var_ptr->dofIndices(),
                                           _vector_tags,
                                           _matrix_tags,
                                           _kappa_var_ptr->scalingFactor());
-#endif
   }
-#else
-  mooseError("residual and jacobian together only supported for global AD indexing");
-#endif
 }
 
 void
@@ -154,56 +127,3 @@ ADKernelScalarBase::computeScalarResidualsForJacobian()
     for (_h = 0; _h < _k_order; _h++)
       _scalar_residuals[_h] += _JxW[_qp] * _coord[_qp] * computeScalarQpResidual();
 }
-
-// void
-// ADKernelScalarBase::computeADScalarJacobian(
-//     const std::vector<std::pair<MooseVariableScalar *, MooseVariableScalar *>> &
-//     coupling_entries)
-// {
-//   computeScalarResidualsForJacobian();
-
-//   auto local_functor =
-//       [&](const std::vector<ADReal> &, const std::vector<dof_id_type> &, const std::set<TagID> &)
-//   {
-//     for (const auto & it : coupling_entries)
-//     {
-//       const MooseVariableScalar & jvariable = *(it.second);
-
-//       if (!jvariable.hasBlocks(_current_elem->subdomain_id()))
-//         continue;
-
-//       // Make sure to get the correct undisplaced/displaced variable
-//       addScalarJacobian(_sys.getScalarVariable(_tid, jvariable.number()));
-//     }
-//   };
-
-//   _assembly.processJacobian(_scalar_residuals,
-//                             dofIndices(),
-//                             _matrix_tags,
-//                             _kappa_var_ptr->scalingFactor(),
-//                             local_functor);
-// }
-
-// void
-// ADKernelScalarBase::addScalarJacobian(const MooseVariableScalar & jvariable)
-// {
-//   unsigned int jvar = jvariable.number();
-
-//   auto ad_offset = Moose::adOffset(jvar, _sys.getMaxVarNDofsPerElem(),
-//   Moose::ElementType::Element);
-
-//   const unsigned int s_order = jvariable.order();
-//   _local_ke.resize(_k_order, s_order);
-
-//   for (_h = 0; _h < _k_order; _h++)
-//     for (_l = 0; _l < s_order; _l++)
-//     {
-// #ifndef MOOSE_SPARSE_AD
-//       mooseAssert(ad_offset + _j < MOOSE_AD_MAX_DOFS_PER_ELEM,
-//                   "Out of bounds access in derivative vector.");
-// #endif
-//       _local_ke(_h, _l) += _scalar_residuals[_h].derivatives()[ad_offset + _l];
-//     }
-
-//   accumulateTaggedLocalMatrix();
-// }

python/TestHarness/testers/PetscJacobianTester.py

@@ -18,7 +18,7 @@ class PetscJacobianTester(RunApp):
     def validParams():
         params = RunApp.validParams()
         params.addParam('ratio_tol', 1e-8, "Relative tolerance to compare the ration against.")
-        params.addParam('difference_tol', 1e-7, "Relative tolerance to compare the difference against.")
+        params.addParam('difference_tol', 1e-6, "Relative tolerance to compare the difference against.")
         params.addParam('state', 'user', "The state for which we want to compare against the "
                                          "finite-differenced Jacobian ('user', 'const_positive', or "
                                          "'const_negative'.")

test/tests/kernels/ad_scalar_kernel_constraint/tests

@@ -1,15 +1,87 @@
 [Tests]
   issues = '#22174'
   design = 'source/kernels/ADScalarLMKernel.md'
-  [./kernel_RJ]
-    type = 'Exodiff'
-    input = 'scalar_constraint_kernel_RJ.i'
-    exodiff = 'scalar_constraint_kernel_RJ_out.e'
-    # This problem only has 4 elements and therefore does not seem to run on > 4 procs.
-    max_parallel = 4
-    requirement = 'The system will produce a harmonic function with a prescribed average value; testing scalar constraint with few other classes.'
-    ad_indexing_type = 'global'
-  [../]
+  [./kernel_dirichlet]
+    requirement = 'The system will produce a harmonic function with a prescribed average value,'
+    [physics_separate]
+      type = 'Exodiff'
+      input = 'scalar_constraint_kernel_RJ.i'
+      exodiff = 'scalar_constraint_kernel_RJ_out.e'
+      # This problem only has 4 elements and therefore does not seem to run on > 4 procs.
+      cli_args = 'Executioner/residual_and_jacobian_together=false'
+      ad_indexing_type = 'global'
+      detail = 'showing the correct results with separate computation of residual and Jacobian'
+    [../]
+    [jacobian_separate]
+      type = 'PetscJacobianTester'
+      run_sim = True
+      input = 'scalar_constraint_kernel_RJ.i'
+      ad_indexing_type = 'global'
+      detail = 'showing the correct results with separate computation of residual and Jacobian'
+      cli_args = 'Executioner/residual_and_jacobian_together=false'
+      prereq = 'kernel_dirichlet/physics_separate'
+    [../]
+    [physics_together]
+      type = 'Exodiff'
+      input = 'scalar_constraint_kernel_RJ.i'
+      exodiff = 'scalar_constraint_kernel_RJ_out.e'
+      # This problem only has 4 elements and therefore does not seem to run on > 4 procs.
+      cli_args = 'Executioner/residual_and_jacobian_together=true'
+      ad_indexing_type = 'global'
+      detail = 'showing the correct results with computation of residual and Jacobian together'
+      prereq = 'kernel_dirichlet/jacobian_separate'
+    [../]
+    [jacobian_together]
+      type = 'PetscJacobianTester'
+      run_sim = True
+      input = 'scalar_constraint_kernel_RJ.i'
+      ad_indexing_type = 'global'
+      detail = 'showing the correct results with computation of residual and Jacobian together'
+      cli_args = 'Executioner/residual_and_jacobian_together=true'
+      prereq = 'kernel_dirichlet/physics_together'
+    [../]
+  []
+
+  [./kernel_neumann]
+    requirement = 'The system shall solve the constrained Neumann problem using the AD Lagrange multiplier approach with LU solver,'
+    [physics_separate]
+      type = 'CSVDiff'
+      input = 'scalar_constraint_together.i'
+      csvdiff = 'scalar_constraint_together_out.csv'
+      # This problem only has 4 elements and therefore does not seem to run on > 4 procs.
+      cli_args = 'Executioner/residual_and_jacobian_together=false'
+      ad_indexing_type = 'global'
+      detail = 'showing the correct results with separate computation of residual and Jacobian'
+    [../]
+    [jacobian_separate]
+      type = 'PetscJacobianTester'
+      run_sim = True
+      input = 'scalar_constraint_together.i'
+      ad_indexing_type = 'global'
+      detail = 'showing the correct results with separate computation of residual and Jacobian'
+      cli_args = 'Executioner/residual_and_jacobian_together=false'
+      prereq = 'kernel_neumann/physics_separate'
+    [../]
+    [physics_together]
+      type = 'CSVDiff'
+      input = 'scalar_constraint_together.i'
+      csvdiff = 'scalar_constraint_together_out.csv'
+      # This problem only has 4 elements and therefore does not seem to run on > 4 procs.
+      cli_args = 'Executioner/residual_and_jacobian_together=true'
+      ad_indexing_type = 'global'
+      detail = 'showing the correct results with computation of residual and Jacobian together'
+      prereq = 'kernel_neumann/jacobian_separate'
+    [../]
+    [jacobian_together]
+      type = 'PetscJacobianTester'
+      run_sim = True
+      input = 'scalar_constraint_together.i'
+      ad_indexing_type = 'global'
+      detail = 'showing the correct results with computation of residual and Jacobian together'
+      cli_args = 'Executioner/residual_and_jacobian_together=true'
+      prereq = 'kernel_neumann/physics_together'
+    [../]
+  []
 
   [./kernel_together]
     type = 'CSVDiff'