Add TEAM4 HForm Mesh Uniform Refinement Test Case

src/auxsolvers/flux_monitor_aux.cpp

@@ -89,6 +89,13 @@ FluxMonitorAux::Init(const hephaestus::GridFunctions & gridfunctions,
   }
 }
 
+void
+FluxMonitorAux::Reset()
+{
+  _times.DeleteAll();
+  _fluxes.DeleteAll();
+}
+
 void
 FluxMonitorAux::Solve(double t)
 {

src/auxsolvers/flux_monitor_aux.hpp

@@ -23,6 +23,8 @@ class FluxMonitorAux : public AuxSolver
   void Init(const hephaestus::GridFunctions & gridfunctions,
             hephaestus::Coefficients & coefficients) override;
 
+  void Reset();
+
   void Update() override {}
 
   void Solve(double t = 0.0) override;

src/sources/scalar_potential_source.cpp

@@ -131,11 +131,8 @@ ScalarPotentialSource::Apply(mfem::ParLinearForm * lf)
   _a0->FormLinearSystem(
       poisson_ess_tdof_list, phi_gf, *_b0, *_diffusion_mat, *_p_tdofs, *_b0_tdofs);
 
-  if (_a0_solver == nullptr)
-  {
-    _a0_solver = std::make_unique<hephaestus::DefaultH1PCGSolver>(_solver_options, *_diffusion_mat);
-  }
   // Solve
+  _a0_solver = std::make_unique<hephaestus::DefaultH1PCGSolver>(_solver_options, *_diffusion_mat);
   _a0_solver->Mult(*_b0_tdofs, *_p_tdofs);
 
   // "undo" the static condensation saving result in grid function dP

test/integration/test_team4_hform.cpp

@@ -43,6 +43,84 @@ class TestTEAM4HForm
     dh_dt(2) = 0.0;
   }
 
+  /// Returns a unique pointer to the constructed problem builder.
+  std::unique_ptr<hephaestus::HFormulation> ConfigureProblem()
+  {
+    // Create Formulation
+    auto problem_builder = std::make_unique<hephaestus::HFormulation>(
+        "electric_resistivity", "electric_conductivity", "magnetic_permeability", "magnetic_field");
+    // Set Mesh
+    mfem::Mesh mesh((std::string(DATA_DIR) + std::string("./team4_symmetrized.g")).c_str(), 1, 1);
+    auto pmesh = std::make_shared<mfem::ParMesh>(MPI_COMM_WORLD, mesh);
+
+    problem_builder->SetMesh(pmesh);
+    problem_builder->AddFESpace("H1", "H1_3D_P1");
+    problem_builder->AddFESpace("HCurl", "ND_3D_P1");
+    problem_builder->AddFESpace("HDiv", "RT_3D_P0");
+    problem_builder->AddFESpace("Scalar_L2", "L2_3D_P0");
+    problem_builder->AddFESpace("Vector_L2", "L2_3D_P0", 3);
+    problem_builder->AddGridFunction("magnetic_field", "HCurl");
+
+    problem_builder->AddGridFunction("dmagnetic_potential_dt", "H1");
+    problem_builder->AddGridFunction("magnetic_flux_density", "HDiv");
+    problem_builder->RegisterMagneticFluxDensityAux("magnetic_flux_density");
+
+    problem_builder->AddGridFunction("current_density", "HDiv");
+    problem_builder->RegisterCurrentDensityAux("current_density");
+
+    problem_builder->AddGridFunction("electric_field", "HCurl");
+    problem_builder->RegisterElectricFieldAux("electric_field");
+
+    hephaestus::Coefficients coefficients = DefineCoefficients();
+    problem_builder->SetCoefficients(coefficients);
+
+    hephaestus::Sources sources = DefineSources();
+    problem_builder->SetSources(sources);
+
+    hephaestus::Outputs outputs = DefineOutputs();
+    problem_builder->SetOutputs(outputs);
+
+    problem_builder->AddBoundaryCondition(
+        "tangential_dhdt_bc",
+        std::make_shared<hephaestus::VectorDirichletBC>(
+            "dmagnetic_field_dt",
+            mfem::Array<int>({1, 2, 5, 6}),
+            coefficients._vectors.Get("surface_tangential_dHdt")));
+    problem_builder->AddBoundaryCondition(
+        "magnetic_potential_bc",
+        std::make_shared<hephaestus::ScalarDirichletBC>(
+            "dmagnetic_potential_dt",
+            mfem::Array<int>({1, 2}),
+            coefficients._scalars.Get("magnetic_potential_time_derivative")));
+
+    auto fluxmonitor = std::make_shared<hephaestus::FluxMonitorAux>("current_density", 3);
+    fluxmonitor->SetPriority(2);
+    problem_builder->AddPostprocessor("FluxMonitor", fluxmonitor);
+
+    hephaestus::InputParameters solver_options;
+    solver_options.SetParam("AbsTolerance", float(1.0e-20));
+    solver_options.SetParam("Tolerance", float(1.0e-20));
+    solver_options.SetParam("MaxIter", (unsigned int)500);
+    problem_builder->SetSolverOptions(solver_options);
+
+    problem_builder->FinalizeProblem();
+
+    return problem_builder;
+  }
+
+  hephaestus::InputParameters DefineExecutionerParameters(hephaestus::TimeDomainProblem & problem)
+  {
+    hephaestus::InputParameters exec_params;
+
+    exec_params.SetParam("TimeStep", float(0.001));
+    exec_params.SetParam("StartTime", float(0.00));
+    exec_params.SetParam("EndTime", float(0.015));
+    exec_params.SetParam("VisualisationSteps", int(1));
+    exec_params.SetParam("Problem", &problem);
+
+    return exec_params;
+  }
+
   hephaestus::Coefficients DefineCoefficients()
   {
     hephaestus::Subdomain brick("brick", 1);
@@ -93,92 +171,105 @@ class TestTEAM4HForm
                      std::make_shared<mfem::ParaViewDataCollection>("Team4ParaView"));
     return outputs;
   }
+
+  /// Sets the peak current and time using the flux monitor.
+  void ExtractPeakCurrentAndTime(hephaestus::TimeDomainProblem & problem,
+                                 double & peak_current,
+                                 double & peak_current_time)
+  {
+    auto flux_monitor = problem._postprocessors.Get<hephaestus::FluxMonitorAux>("FluxMonitor");
+
+    const mfem::real_t min_flux = flux_monitor->_fluxes.Min();
+
+    peak_current = -2.0 * min_flux;
+    peak_current_time = 0.0;
+
+    for (int i = 0; i < flux_monitor->_times.Size(); i++)
+    {
+      const mfem::real_t flux = flux_monitor->_fluxes[i];
+      const mfem::real_t flux_time = flux_monitor->_times[i];
+
+      if (flux == min_flux)
+      {
+        peak_current_time = flux_time;
+      }
+    }
+  }
+
+  /// Checks peak current and time are within expected range.
+  void VerifyPeakCurrentAndTimeWithinTolerances(double & peak_current, double & peak_current_time)
+  {
+    const double min_peak_current = 3.2e3;
+    const double max_peak_current = 3.6e3;
+
+    const double min_peak_current_time = 0.010;
+    const double max_peak_current_time = 0.012;
+
+    REQUIRE(peak_current > min_peak_current);
+    REQUIRE(peak_current < max_peak_current);
+
+    REQUIRE(peak_current_time > min_peak_current_time);
+    REQUIRE(peak_current_time < max_peak_current_time);
+  }
+
+  /// Calls Reset method of FluxMonitorAux.
+  void ResetFluxMonitor(const hephaestus::TimeDomainProblem & problem)
+  {
+    auto fluxmonitor = problem._postprocessors.Get<hephaestus::FluxMonitorAux>("FluxMonitor");
+
+    fluxmonitor->Reset();
+  }
 };
 
 TEST_CASE_METHOD(TestTEAM4HForm, "TestTEAM4HForm", "[CheckRun]")
 {
-  // Create Formulation
-  auto problem_builder = std::make_unique<hephaestus::HFormulation>(
-      "electric_resistivity", "electric_conductivity", "magnetic_permeability", "magnetic_field");
-  // Set Mesh
-  mfem::Mesh mesh((std::string(DATA_DIR) + std::string("./team4_symmetrized.g")).c_str(), 1, 1);
-  auto pmesh = std::make_shared<mfem::ParMesh>(MPI_COMM_WORLD, mesh);
-
-  problem_builder->SetMesh(pmesh);
-  problem_builder->AddFESpace("H1", "H1_3D_P1");
-  problem_builder->AddFESpace("HCurl", "ND_3D_P1");
-  problem_builder->AddFESpace("HDiv", "RT_3D_P0");
-  problem_builder->AddFESpace("Scalar_L2", "L2_3D_P0");
-  problem_builder->AddFESpace("Vector_L2", "L2_3D_P0", 3);
-  problem_builder->AddGridFunction("magnetic_field", "HCurl");
-
-  problem_builder->AddGridFunction("dmagnetic_potential_dt", "H1");
-  problem_builder->AddGridFunction("magnetic_flux_density", "HDiv");
-  problem_builder->RegisterMagneticFluxDensityAux("magnetic_flux_density");
-
-  problem_builder->AddGridFunction("current_density", "HDiv");
-  problem_builder->RegisterCurrentDensityAux("current_density");
-
-  problem_builder->AddGridFunction("electric_field", "HCurl");
-  problem_builder->RegisterElectricFieldAux("electric_field");
-
-  hephaestus::Coefficients coefficients = DefineCoefficients();
-  problem_builder->SetCoefficients(coefficients);
-
-  hephaestus::Sources sources = DefineSources();
-  problem_builder->SetSources(sources);
-
-  hephaestus::Outputs outputs = DefineOutputs();
-  problem_builder->SetOutputs(outputs);
-
-  problem_builder->AddBoundaryCondition("tangential_dhdt_bc",
-                                        std::make_shared<hephaestus::VectorDirichletBC>(
-                                            "dmagnetic_field_dt",
-                                            mfem::Array<int>({1, 2, 5, 6}),
-                                            coefficients._vectors.Get("surface_tangential_dHdt")));
-  problem_builder->AddBoundaryCondition(
-      "magnetic_potential_bc",
-      std::make_shared<hephaestus::ScalarDirichletBC>(
-          "dmagnetic_potential_dt",
-          mfem::Array<int>({1, 2}),
-          coefficients._scalars.Get("magnetic_potential_time_derivative")));
-
-  auto fluxmonitor = std::make_shared<hephaestus::FluxMonitorAux>("current_density", 3);
-  fluxmonitor->SetPriority(2);
-  problem_builder->AddPostprocessor("FluxMonitor", fluxmonitor);
-
-  hephaestus::InputParameters solver_options;
-  solver_options.SetParam("AbsTolerance", float(1.0e-20));
-  solver_options.SetParam("Tolerance", float(1.0e-20));
-  solver_options.SetParam("MaxIter", (unsigned int)500);
-  problem_builder->SetSolverOptions(solver_options);
-
-  problem_builder->FinalizeProblem();
-
+  auto problem_builder = ConfigureProblem();
   auto problem = problem_builder->ReturnProblem();
-  hephaestus::InputParameters exec_params;
-  exec_params.SetParam("TimeStep", float(0.001));
-  exec_params.SetParam("StartTime", float(0.00));
-  exec_params.SetParam("EndTime", float(0.015));
-  exec_params.SetParam("VisualisationSteps", int(1));
-  exec_params.SetParam("Problem", static_cast<hephaestus::TimeDomainProblem *>(problem.get()));
 
-  auto executioner = std::make_unique<hephaestus::TransientExecutioner>(exec_params);
+  hephaestus::InputParameters exec_params = DefineExecutionerParameters(*problem);
+
+  hephaestus::TransientExecutioner executioner(exec_params);
+  executioner.Execute();
+
+  double peak_current, peak_current_time;
+  ExtractPeakCurrentAndTime(*problem, peak_current, peak_current_time);
+
+  VerifyPeakCurrentAndTimeWithinTolerances(peak_current, peak_current_time);
+}
 
-  executioner->Execute();
+/// Test case for checking "Updates" work as expected following a mesh refinement.
+TEST_CASE_METHOD(TestTEAM4HForm, "TestTEAM4HFormMeshUpdates", "[CheckRun][!benchmark]")
+{
+  auto problem_builder = ConfigureProblem();
+  auto problem = problem_builder->ReturnProblem();
 
-  double peak_current = -2 * fluxmonitor->_fluxes.Min();
-  double peak_current_time;
-  for (std::size_t i = 0; i < fluxmonitor->_times.Size(); ++i)
+  hephaestus::InputParameters exec_params = DefineExecutionerParameters(*problem);
+  hephaestus::logger.set_level(spdlog::level::info);
+
+  const int imax_refinement = 2;
+  for (int irefinement = 0; irefinement < imax_refinement; irefinement++)
   {
-    if (fluxmonitor->_fluxes[i] == fluxmonitor->_fluxes.Min())
+    // Refine and reset flux monitor on subsequent runs.
+    if (irefinement != 0)
     {
-      peak_current_time = fluxmonitor->_times[i];
+      problem->_pmesh->UniformRefinement();
+      problem->Update();
+
+      ResetFluxMonitor(*problem);
     }
-  }
 
-  REQUIRE(peak_current > 3.2e3);
-  REQUIRE(peak_current < 3.6e3);
-  REQUIRE(peak_current_time < 0.012);
-  REQUIRE(peak_current_time > 0.01);
+    // Create a new executioner.
+    hephaestus::TransientExecutioner executioner(exec_params);
+    executioner.Execute();
+
+    double peak_current, peak_current_time;
+    ExtractPeakCurrentAndTime(*problem, peak_current, peak_current_time);
+
+    hephaestus::logger.info("Refinement Level: {}, Peak current: {}, Time: {} seconds",
+                            irefinement,
+                            peak_current,
+                            peak_current_time);
+
+    VerifyPeakCurrentAndTimeWithinTolerances(peak_current, peak_current_time);
+  }
 }