Skip to content

Add TEAM4 HForm Mesh Uniform Refinement Test Case #113

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
EdwardPalmer99 wants to merge 14 commits into
base: EdwardPalmer99/mesh-update-feature-branch
Choose a base branch
from
Open

Add TEAM4 HForm Mesh Uniform Refinement Test Case #113

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
9a7b843
Added Reset method to FluxMonitorAux.
May 20, 2024
135763d
Added additional mesh update test case to test_team4_Hform which runs…
May 20, 2024
f53313b
Fixes issue with test4 where the serial mesh was mistakenly refined r…
May 24, 2024
3166048
Rebuild a0 solver in ScalarPotentialSource when Apply called.
May 24, 2024
af4d36f
Merge branch 'EdwardPalmer99/mesh-update-feature/rebuild-scalar-poten…
May 24, 2024
53ffaa1
Reduces max refinement level due to time taken to execute problem.
May 24, 2024
b73a233
Adds "[!benchmark]" tag to TestTeam4HFormMeshUpdates to ensure it is …
May 30, 2024
195e67c
Adds "ConfigureProblem" method to TestTEAM4HForm to reduce code dupli…
May 30, 2024
f344c4d
Adds DefineExecutionerParameters method.
May 30, 2024
617ee80
Adds ExtractPeakCurrentAndTime method.
May 30, 2024
7112196
Adds "VerifyPeakCurrentAndTimeWithinTolerances".
May 30, 2024
440421d
Updates comments in TestTEAM4HForm.
May 30, 2024
670cce9
Moves pmesh refinement to top of for loop.
May 30, 2024
bc3a110
Removes use of unique pointer for TransientExecutioner in TestTEAM4HF…
May 30, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions src/auxsolvers/flux_monitor_aux.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,13 @@ FluxMonitorAux::Init(const hephaestus::GridFunctions & gridfunctions,
}
}

void
FluxMonitorAux::Reset()
{
_times.DeleteAll();
_fluxes.DeleteAll();
}

void
FluxMonitorAux::Solve(double t)
{
Expand Down
2 changes: 2 additions & 0 deletions src/auxsolvers/flux_monitor_aux.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
5 changes: 1 addition & 4 deletions src/sources/scalar_potential_source.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Copy link
Contributor

@alexanderianblair alexanderianblair Jul 12, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Duplicate of #112 ?

_a0_solver->Mult(*_b0_tdofs, *_p_tdofs);

// "undo" the static condensation saving result in grid function dP
Expand Down
243 changes: 167 additions & 76 deletions test/integration/test_team4_hform.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
}
}