Merge pull request #3421 from roystgnr/still_more_coverage

Still more test coverage

Author: roystgnr

examples/adaptivity/adaptivity_ex4/adaptivity_ex4.C

@@ -58,6 +58,8 @@
 #include "libmesh/perf_log.h"
 #include "libmesh/string_to_enum.h"
 #include "libmesh/enum_solver_package.h"
+#include "libmesh/dirichlet_boundaries.h"
+#include "libmesh/wrapped_function.h"
 
 // Bring in everything from the libMesh namespace
 using namespace libMesh;
@@ -71,6 +73,9 @@ using namespace libMesh;
 void assemble_biharmonic(EquationSystems & es,
                          const std::string & system_name);
 
+// Global parameter, to be accessible to the above function
+bool penalty_dirichlet = true;
+
 
 // Prototypes for calculation of the exact solution.  Necessary
 // for setting boundary conditions.
@@ -189,6 +194,7 @@ int main(int argc, char ** argv)
   const Real coarsen_percentage = input_file("coarsen_percentage", 0.5);
   const unsigned int dim = input_file("dimension", 2);
   const unsigned int max_linear_iterations = input_file("max_linear_iterations", 10000);
+  penalty_dirichlet = input_file("penalty_dirichlet", penalty_dirichlet);
 
   // Skip higher-dimensional examples on a lower-dimensional libMesh build
   libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");
@@ -217,10 +223,12 @@ int main(int argc, char ** argv)
   else
     output_file += "clough_";
 
+  output_file += approx_order_string;
+
   if (uniform_refine == 0)
-    output_file += "adaptive";
+    output_file += "_adaptive";
   else
-    output_file += "uniform";
+    output_file += "_uniform";
 
 #ifdef LIBMESH_HAVE_EXODUS_API
   // If we have Exodus, use the same base output filename
@@ -303,6 +311,26 @@ int main(int argc, char ** argv)
   // function.
   system.attach_assemble_function (assemble_biharmonic);
 
+  // Give the system Dirichlet boundary conditions, if we want to
+  // enforce those strongly.
+  if (!penalty_dirichlet)
+    {
+      const std::vector<unsigned int>
+        all_vars(1, system.variable_number("u"));
+      std::set<boundary_id_type> all_bdys;
+      for (unsigned int i=0; i != dim; ++i)
+        {
+          all_bdys.insert(2*i);
+          all_bdys.insert(2*i+1);
+        }
+
+      WrappedFunction<Number> exact_val(system, exact_solution);
+      WrappedFunction<Gradient> exact_grad(system, *exact_derivative);
+      DirichletBoundary exact_bc(all_bdys, all_vars, exact_val,
+                                 exact_grad);
+      system.get_dof_map().add_dirichlet_boundary(exact_bc);
+    }
+
   // Initialize the data structures for the equation system.
   equation_systems.init();
 
@@ -377,6 +405,20 @@ int main(int argc, char ** argv)
           << exact_sol.h1_error("Biharmonic", "u") << " "
           << exact_sol.h2_error("Biharmonic", "u") << std::endl;
 
+      // Assert that we're not outright diverging here, for regression
+      // testing.  These bounds just come from eyeballing some graphs
+      // and tacking on a huge tolerance.
+      const dof_id_type n_dofs = system.n_dofs();
+
+      if (exact_sol.l2_error("Biharmonic", "u") > 200*std::pow(n_dofs,-4./3.))
+        libmesh_error_msg("Unacceptably high L2-error!");
+
+      if (exact_sol.h1_error("Biharmonic", "u") > 200./n_dofs)
+        libmesh_error_msg("Unacceptably high H1-error!");
+
+      if (exact_sol.h2_error("Biharmonic", "u") > 400*std::pow(n_dofs,-2./3.))
+        libmesh_error_msg("Unacceptably high H2-error!");
+
       // Possibly refine the mesh
       if (r_step+1 != max_r_steps)
         {
@@ -836,12 +878,14 @@ void assemble_biharmonic(EquationSystems & es,
 
 
       // At this point the interior element integration has
-      // been completed.  However, we have not yet addressed
-      // boundary conditions.  For this example we will only
-      // consider simple Dirichlet boundary conditions imposed
+      // been completed.  However, we may have not yet addressed
+      // boundary conditions, if we didn't add a DirichletBoundary to
+      // our system earlier.  In that case, we will demonstrate the
+      // imposition of simple Dirichlet boundary conditions imposed
       // via the penalty method.  Note that this is a fourth-order
       // problem: Dirichlet boundary conditions include *both*
       // boundary values and boundary normal fluxes.
+      if (penalty_dirichlet)
       {
         // Start logging the boundary condition computation.  We use a
         // macro to log everything in this scope.

examples/adaptivity/adaptivity_ex4/run.sh

@@ -7,12 +7,15 @@ source $LIBMESH_DIR/examples/run_common.sh
 example_name=adaptivity_ex4
 
 run_example "$example_name" approx_type=HERMITE
+run_example "$example_name" approx_type=HERMITE penalty_dirichlet=false
 
 # HERMITE elements can p refine in 1D
 run_example "$example_name" dimension=1 approx_type=HERMITE approx_order=FOURTH
+run_example "$example_name" dimension=1 approx_type=HERMITE approx_order=FOURTH penalty_dirichlet=false
 
 # CLOUGH elements don't currently support threads
 run_example_no_extra_options "$example_name" approx_type=CLOUGH $LIBMESH_OPTIONS --n_threads=1
+run_example_no_extra_options "$example_name" approx_type=CLOUGH $LIBMESH_OPTIONS --n_threads=1 penalty_dirichlet=false
 #run_example "$example_name" approx_type=SECOND  # Broken?
 
 # Examples to use for benchmarking

examples/adjoints/adjoints_ex2/adjoints_ex2.C

@@ -509,12 +509,10 @@ int main (int argc, char ** argv)
 
         write_output(equation_systems, a_step, "primal", param);
 
-        NumericVector<Number> & primal_solution = *system.solution;
+        system.assemble_qoi_sides = true;
 
         SensitivityData sensitivities(qois, system, system.get_parameter_vector());
 
-        system.assemble_qoi_sides = true;
-
         // Here we solve the adjoint problem inside the adjoint_qoi_parameter_sensitivity
         // function, so we have to set the adjoint_already_solved boolean to false
         system.set_adjoint_already_solved(false);
@@ -563,8 +561,45 @@ int main (int argc, char ** argv)
         libmesh_assert_less(std::abs((sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact), 2.e-4);
         libmesh_assert_less(std::abs((sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact), 2.e-4);
 
-        NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
+        // Let's do a forward sensitivity solve too, unless we're
+        // told to skip it for backwards compatibility with old
+        // performance benchmarks.
+        const bool forward_sensitivity = infile("--forward_sensitivity", true);
+
+        if (forward_sensitivity)
+          {
+            // This will require two linear solves (one per parameter)
+            // rather than the adjoint sensitivity's one, but it's useful
+            // for regression testing.
+            SensitivityData forward_sensitivities(qois, system, system.get_parameter_vector());
+            system.forward_qoi_parameter_sensitivity(qois, system.get_parameter_vector(), forward_sensitivities);
+
+            libmesh_assert_less(std::abs((forward_sensitivities[0][0] - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact), 2.e-4);
+            libmesh_assert_less(std::abs((forward_sensitivities[0][1] - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact), 2.e-4);
+
+            // These should be the same linearization, just calculated
+            // different ways with different roundoff error
+            libmesh_assert_less
+              (std::abs((forward_sensitivities[0][0] - sensitivity_QoI_0_0_computed)/sensitivity_QoI_0_0_computed), TOLERANCE);
+            libmesh_assert_less
+              (std::abs((forward_sensitivities[0][1] - sensitivity_QoI_0_1_computed)/sensitivity_QoI_0_1_computed), TOLERANCE);
+
+            libMesh::out << "The error in forward calculation of sensitivity QoI_0_0 is "
+                         << std::setprecision(17)
+                         << std::abs(forward_sensitivities[0][0] - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact
+                         << std::endl;
+
+            libMesh::out << "The error in forward calculation of sensitivity QoI_0_1 is "
+                         << std::setprecision(17)
+                         << std::abs(forward_sensitivities[0][1] - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact
+                         << std::endl
+                         << std::endl;
+
+
+          }
 
+        NumericVector<Number> & primal_solution = *system.solution;
+        NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
         primal_solution.swap(dual_solution_0);
         write_output(equation_systems, a_step, "adjoint_0", param);
 

examples/adjoints/adjoints_ex2/run.sh

@@ -11,4 +11,4 @@ example_dir=examples/adjoints/$example_name
 run_example "$example_name"
 
 # Benchmark parameters
-benchmark_example 1 "$example_name" coarserefinements=4
+benchmark_example 1 "$example_name" coarserefinements=4 --forward_sensitivity=false

examples/fem_system/fem_system_ex1/run.sh

@@ -9,7 +9,9 @@ example_name=fem_system_ex1
 example_dir=examples/fem_system/$example_name
 
 # First, run with standard options from input files.
-run_example "$example_name"
+# We'll throw Jacobian verification in too, just to get some test
+# coverage there.
+run_example "$example_name" verify_analytic_jacobians=1e-6
 
 # Next, lets run with pure fieldsplit without gmg, if we have PETSc
 

include/partitioning/centroid_partitioner.h

@@ -82,6 +82,8 @@ class CentroidPartitioner : public Partitioner
   CentroidPartitioner & operator= (CentroidPartitioner &&) = default;
   virtual ~CentroidPartitioner() = default;
 
+  virtual PartitionerType type () const override;
+
   /**
    * \returns A copy of this partitioner wrapped in a smart pointer.
    */

include/partitioning/hilbert_sfc_partitioner.h

@@ -23,6 +23,8 @@
 // Local Includes
 #include "libmesh/sfc_partitioner.h"
 
+#include "libmesh/enum_partitioner_type.h"
+
 namespace libMesh
 {
 
@@ -56,6 +58,11 @@ class HilbertSFCPartitioner : public SFCPartitioner
   HilbertSFCPartitioner & operator= (HilbertSFCPartitioner &&) = default;
   virtual ~HilbertSFCPartitioner() = default;
 
+  virtual PartitionerType type () const override
+  {
+    return HILBERT_SFC_PARTITIONER;
+  };
+
   /**
    * \returns A copy of this partitioner wrapped in a smart pointer.
    */

include/partitioning/linear_partitioner.h

@@ -54,6 +54,8 @@ class LinearPartitioner : public Partitioner
   LinearPartitioner & operator= (LinearPartitioner &&) = default;
   virtual ~LinearPartitioner() = default;
 
+  virtual PartitionerType type () const override;
+
   /**
    * \returns A copy of this partitioner wrapped in a smart pointer.
    */

include/partitioning/mapped_subdomain_partitioner.h

@@ -56,6 +56,8 @@ class MappedSubdomainPartitioner : public Partitioner
   MappedSubdomainPartitioner & operator= (MappedSubdomainPartitioner &&) = default;
   virtual ~MappedSubdomainPartitioner() = default;
 
+  virtual PartitionerType type () const override;
+
   /**
    * \returns A copy of this partitioner wrapped in a smart pointer.
    */

include/partitioning/metis_partitioner.h

@@ -52,6 +52,8 @@ class MetisPartitioner : public Partitioner
   MetisPartitioner & operator= (MetisPartitioner &&) = default;
   virtual ~MetisPartitioner() = default;
 
+  virtual PartitionerType type () const override;
+
   /**
    * \returns A copy of this partitioner wrapped in a smart pointer.
    */