Merge pull request #3519 from roystgnr/residual_constrain_on_elements

Apply constraints via element nonlinear residuals

Author: roystgnr

examples/adjoints/adjoints_ex1/adjoints_ex1.C

@@ -67,6 +67,7 @@
 #include "libmesh/newton_solver.h"
 #include "libmesh/numeric_vector.h"
 #include "libmesh/partitioner.h"
+#include "libmesh/petsc_diff_solver.h"
 #include "libmesh/steady_solver.h"
 #include "libmesh/system_norm.h"
 #include "libmesh/enum_solver_package.h"
@@ -184,6 +185,16 @@ void set_system_parameters(LaplaceSystem & system,
   system.time_solver = std::make_unique<SteadySolver>(system);
 
   // Nonlinear solver options
+  if (param.use_petsc_snes)
+  {
+#ifdef LIBMESH_HAVE_PETSC
+    PetscDiffSolver *solver = new PetscDiffSolver(system);
+    system.time_solver->diff_solver().reset(solver);
+#else
+    libmesh_error_msg("This example requires libMesh to be compiled with PETSc support.");
+#endif
+  }
+  else
   {
     NewtonSolver * solver = new NewtonSolver(system);
     system.time_solver->diff_solver() = std::unique_ptr<DiffSolver>(solver);
@@ -201,6 +212,7 @@ void set_system_parameters(LaplaceSystem & system,
         solver->continue_after_max_iterations = true;
         solver->continue_after_backtrack_failure = true;
       }
+    system.set_constrain_in_solver(param.constrain_in_solver);
 
     // And the linear solver options
     solver->max_linear_iterations       = param.max_linear_iterations;

examples/adjoints/adjoints_ex1/femparameters.C

@@ -81,6 +81,7 @@ FEMParameters::FEMParameters(const Parallel::Communicator & comm_in) :
   print_element_residuals(false),
   print_element_jacobians(false),
 
+  constrain_in_solver(true),
   use_petsc_snes(false),
   time_solver_quiet(true), solver_quiet(true), solver_verbose(false),
   reuse_preconditioner(true),
@@ -686,6 +687,7 @@ void FEMParameters::read(GetPot & input,
   GETPOT_INPUT(print_element_jacobians);
 
 
+  GETPOT_INPUT(constrain_in_solver);
   GETPOT_INPUT(use_petsc_snes);
   GETPOT_INPUT(time_solver_quiet);
   GETPOT_INPUT(solver_quiet);

examples/adjoints/adjoints_ex1/femparameters.h

@@ -125,6 +125,7 @@ class FEMParameters : public libMesh::ParallelObject
 
   //   Solver options
 
+  bool constrain_in_solver;
   bool use_petsc_snes;
   bool time_solver_quiet, solver_quiet, solver_verbose,
     reuse_preconditioner, require_residual_reduction;

examples/adjoints/adjoints_ex2/adjoints_ex2.C

@@ -69,6 +69,7 @@
 #include "libmesh/mesh_refinement.h"
 #include "libmesh/newton_solver.h"
 #include "libmesh/numeric_vector.h"
+#include "libmesh/petsc_diff_solver.h"
 #include "libmesh/steady_solver.h"
 #include "libmesh/system_norm.h"
 #include "libmesh/enum_solver_package.h"
@@ -187,6 +188,16 @@ void set_system_parameters(LaplaceSystem & system, FEMParameters & param)
   system.time_solver = std::make_unique<SteadySolver>(system);
 
   // Nonlinear solver options
+  if (param.use_petsc_snes)
+  {
+#ifdef LIBMESH_HAVE_PETSC
+    PetscDiffSolver *solver = new PetscDiffSolver(system);
+    system.time_solver->diff_solver().reset(solver);
+#else
+    libmesh_error_msg("This example requires libMesh to be compiled with PETSc support.");
+#endif
+  }
+  else
   {
     NewtonSolver * solver = new NewtonSolver(system);
     system.time_solver->diff_solver() = std::unique_ptr<DiffSolver>(solver);
@@ -203,6 +214,7 @@ void set_system_parameters(LaplaceSystem & system, FEMParameters & param)
         solver->continue_after_max_iterations = true;
         solver->continue_after_backtrack_failure = true;
       }
+    system.set_constrain_in_solver(param.constrain_in_solver);
 
     // And the linear solver options
     solver->max_linear_iterations    = param.max_linear_iterations;

examples/adjoints/adjoints_ex2/femparameters.h

@@ -125,6 +125,7 @@ class FEMParameters : public libMesh::ParallelObject
 
   //   Solver options
 
+  bool constrain_in_solver;
   bool use_petsc_snes;
   bool time_solver_quiet, solver_quiet, solver_verbose,
     reuse_preconditioner, require_residual_reduction;

examples/adjoints/adjoints_ex3/adjoints_ex3.C

@@ -579,6 +579,7 @@ void set_system_parameters(FEMSystem & system,
           solver->continue_after_max_iterations = true;
           solver->continue_after_backtrack_failure = true;
         }
+      system.set_constrain_in_solver(param.constrain_in_solver);
 
       // And the linear solver options
       solver->max_linear_iterations       = param.max_linear_iterations;

examples/adjoints/adjoints_ex3/femparameters.h

@@ -125,6 +125,7 @@ class FEMParameters : public libMesh::ParallelObject
 
   //   Solver options
 
+  bool constrain_in_solver;
   bool use_petsc_snes;
   bool time_solver_quiet, solver_quiet, solver_verbose,
     reuse_preconditioner, require_residual_reduction;

examples/adjoints/adjoints_ex4/adjoints_ex4.C

@@ -65,6 +65,7 @@
 #include "libmesh/mesh_refinement.h"
 #include "libmesh/newton_solver.h"
 #include "libmesh/numeric_vector.h"
+#include "libmesh/petsc_diff_solver.h"
 #include "libmesh/steady_solver.h"
 #include "libmesh/system_norm.h"
 #include "libmesh/enum_solver_package.h"
@@ -177,6 +178,16 @@ void set_system_parameters(LaplaceSystem & system,
   system.time_solver = std::make_unique<SteadySolver>(system);
 
   // Nonlinear solver options
+  if (param.use_petsc_snes)
+  {
+#ifdef LIBMESH_HAVE_PETSC
+    PetscDiffSolver *solver = new PetscDiffSolver(system);
+    system.time_solver->diff_solver().reset(solver);
+#else
+    libmesh_error_msg("This example requires libMesh to be compiled with PETSc support.");
+#endif
+  }
+  else
   {
     NewtonSolver * solver = new NewtonSolver(system);
     system.time_solver->diff_solver() = std::unique_ptr<DiffSolver>(solver);
@@ -193,6 +204,7 @@ void set_system_parameters(LaplaceSystem & system,
         solver->continue_after_max_iterations = true;
         solver->continue_after_backtrack_failure = true;
       }
+    system.set_constrain_in_solver(param.constrain_in_solver);
 
     // And the linear solver options
     solver->max_linear_iterations       = param.max_linear_iterations;

examples/adjoints/adjoints_ex4/femparameters.h

@@ -125,6 +125,7 @@ class FEMParameters : public libMesh::ParallelObject
 
   //   Solver options
 
+  bool constrain_in_solver;
   bool use_petsc_snes;
   bool time_solver_quiet, solver_quiet, solver_verbose,
     reuse_preconditioner, require_residual_reduction;

examples/adjoints/adjoints_ex5/adjoints_ex5.C

@@ -98,6 +98,7 @@
 #include "libmesh/euler_solver.h"
 #include "libmesh/euler2_solver.h"
 #include "libmesh/newton_solver.h"
+#include "libmesh/petsc_diff_solver.h"
 #include "libmesh/twostep_time_solver.h"
 
 #include "libmesh/getpot.h"
@@ -338,6 +339,7 @@ void set_system_parameters(HeatSystem &system, FEMParameters &param)
           solver->continue_after_max_iterations = true;
           solver->continue_after_backtrack_failure = true;
         }
+      system.set_constrain_in_solver(param.constrain_in_solver);
 
       // And the linear solver options
       solver->max_linear_iterations       = param.max_linear_iterations;