diff --git a/data/hollow_sphere_vac.e b/data/hollow_sphere_vac.e new file mode 100644 index 000000000..43777c8dc Binary files /dev/null and b/data/hollow_sphere_vac.e differ diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 7f1b42806..404c9f868 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -4,4 +4,5 @@ add_subdirectory(team7) add_subdirectory(complex_team7) add_subdirectory(closed_coil) add_subdirectory(open_coil) -add_subdirectory(magnetostatic) \ No newline at end of file +add_subdirectory(magnetostatic) +add_subdirectory(hollow_sphere_ren1994) \ No newline at end of file diff --git a/examples/hollow_sphere_ren1994/CMakeLists.txt b/examples/hollow_sphere_ren1994/CMakeLists.txt new file mode 100644 index 000000000..1e43e980e --- /dev/null +++ b/examples/hollow_sphere_ren1994/CMakeLists.txt @@ -0,0 +1,24 @@ +file(GLOB_RECURSE example_src "*.cpp") + +file(GLOB_RECURSE test_src_files "${PROJECT_SOURCE_DIR}/src/*.h" "${PROJECT_SOURCE_DIR}/src/*.hpp" "${PROJECT_SOURCE_DIR}/src/*.cpp") + +set (${PROJECT_NAME}_INCLUDE_DIRS "") +foreach (_srcFile ${test_src_files}) + get_filename_component(_dir ${_srcFile} PATH) + list (APPEND ${PROJECT_NAME}_INCLUDE_DIRS ${_dir}) +endforeach() +list (REMOVE_DUPLICATES ${PROJECT_NAME}_INCLUDE_DIRS) + + +add_executable(hollow_sphere_ren1994 ${example_src}) +add_compile_options(hollow_sphere_ren1994 ${BUILD_TYPE_COMPILER_FLAGS}) +target_include_directories(hollow_sphere_ren1994 PUBLIC ${MFEM_COMMON_INCLUDES} ${MFEM_INCLUDE_DIRS}) +target_include_directories(hollow_sphere_ren1994 PUBLIC ${${PROJECT_NAME}_INCLUDE_DIRS}) +target_include_directories(hollow_sphere_ren1994 PUBLIC ${PROJECT_SOURCE_DIR}/data/) + +set_property(TARGET hollow_sphere_ren1994 PROPERTY CXX_STANDARD 17) + +target_link_libraries(hollow_sphere_ren1994 spdlog::spdlog) +target_link_libraries(hollow_sphere_ren1994 ${GTEST_LIBRARY} ${GTEST_MAIN_LIBRARY} pthread) +target_link_libraries(hollow_sphere_ren1994 ${CMAKE_LIBRARY_OUTPUT_DIRECTORY}/lib${PROJECT_NAME}.so) +target_link_libraries(hollow_sphere_ren1994 ${MFEM_LIBRARIES} ${MFEM_COMMON_LIBRARY} -lrt) diff --git a/examples/hollow_sphere_ren1994/Main.cpp b/examples/hollow_sphere_ren1994/Main.cpp new file mode 100644 index 000000000..ecc88163a --- /dev/null +++ b/examples/hollow_sphere_ren1994/Main.cpp @@ -0,0 +1,208 @@ +#include "hephaestus.hpp" + +const char * DATA_DIR = "../data/"; +const char * MESH_NAME = "hollow_sphere_vac.e"; + +hephaestus::Coefficients +defineCoefficients() +{ + double air_permeability = M_PI * 4.0e-7; + double solid_permeability = air_permeability * 500.0; + + double solid_conductivity = 1.0; + double air_conductivity = 1.0; + + hephaestus::Subdomain vacuum_region("vacuum_region", 107); + vacuum_region._scalar_coefficients.Register( + "electrical_conductivity", std::make_shared(air_conductivity)); + vacuum_region._scalar_coefficients.Register( + "magnetic_permeability", std::make_shared(air_permeability)); + + hephaestus::Subdomain sphere("sphere", 100); + sphere._scalar_coefficients.Register( + "electrical_conductivity", std::make_shared(solid_conductivity)); + sphere._scalar_coefficients.Register( + "magnetic_permeability", std::make_shared(solid_permeability)); + + hephaestus::Subdomain coil_0("coil_0", 103); + coil_0._scalar_coefficients.Register( + "electrical_conductivity", std::make_shared(air_conductivity)); + coil_0._scalar_coefficients.Register( + "magnetic_permeability", std::make_shared(air_permeability)); + + hephaestus::Subdomain coil_1("coil_1", 104); + coil_1._scalar_coefficients.Register( + "electrical_conductivity", std::make_shared(air_conductivity)); + coil_1._scalar_coefficients.Register( + "magnetic_permeability", std::make_shared(air_permeability)); + + hephaestus::Subdomain coil_2("coil_2", 105); + coil_2._scalar_coefficients.Register( + "electrical_conductivity", std::make_shared(air_conductivity)); + coil_2._scalar_coefficients.Register( + "magnetic_permeability", std::make_shared(air_permeability)); + + hephaestus::Subdomain coil_3("coil_3", 106); + coil_3._scalar_coefficients.Register( + "electrical_conductivity", std::make_shared(air_conductivity)); + coil_3._scalar_coefficients.Register( + "magnetic_permeability", std::make_shared(air_permeability)); + + hephaestus::Coefficients coefficients( + std::vector({vacuum_region, sphere, coil_0, coil_1, coil_2, coil_3})); + // coefficients._scalars.Register("frequency", + // std::make_shared(200.0)); + // coefficients._scalars.Register("dielectric_permittivity", + // std::make_shared(8.854e-12)); + + coefficients._scalars.Register("I", std::make_shared(20000)); + + return coefficients; +} + +hephaestus::Sources +defineSources() +{ + hephaestus::InputParameters div_free_source_params; + + // This vector of subdomains will form the coil that we pass to + // ClosedCoilSolver + int order = 1; + int electrode_attr = 110; // coil face + std::string coil_attr = "103 104 105 106"; + mfem::Array coil_domains; + std::stringstream ss(coil_attr); + int att; + + while (ss >> att) + { + coil_domains.Append(att); + std::cout << "coil domains appending " << att << std::endl; + } + + hephaestus::Sources sources; + sources.Register("source", + std::make_shared("source_grad_phi", + "HCurl", + "H1", + "I", + "electrical_conductivity", + coil_domains, + electrode_attr, + true)); + return sources; +} + +hephaestus::Outputs +defineOutputs() +{ + hephaestus::Outputs outputs; + outputs.Register("ParaViewDataCollection", + std::make_shared("HollowSphereParaView")); + return outputs; +} + +int +main(int argc, char * argv[]) +{ + int ref_level(0); + mfem::OptionsParser args(argc, argv); + args.AddOption( + &DATA_DIR, "-dataDir", "--data_directory", "Directory storing input data for tests."); + args.AddOption(&MESH_NAME, "-meshName", "--mesh_name", "File name."); + args.AddOption( + &ref_level, "-ref", "--refinement_level", "number of uniform refinement iterations."); + args.Parse(); + MPI_Init(&argc, &argv); + + hephaestus::logger.set_level(spdlog::level::info); + + // Create Formulation + auto problem_builder = std::make_unique( + "magnetic_reluctivity", "magnetic_permeability", "magnetic_vector_potential"); + + // Set Mesh + mfem::Mesh mesh((std::string(DATA_DIR) + std::string(MESH_NAME)).c_str(), 1, 1); + auto pmesh = std::make_shared(MPI_COMM_WORLD, mesh); + + for (int l = 0; l < ref_level; ++l) + pmesh->UniformRefinement(); + + problem_builder->SetMesh(pmesh); + problem_builder->AddFESpace(std::string("H1"), std::string("H1_3D_P1")); + problem_builder->AddFESpace(std::string("Vector_H1"), std::string("H1_3D_P1"), 3); + problem_builder->AddFESpace(std::string("HCurl"), std::string("ND_3D_P1")); + problem_builder->AddFESpace(std::string("HDiv"), std::string("RT_3D_P0")); + problem_builder->AddFESpace(std::string("Scalar_L2"), std::string("L2_3D_P0")); + problem_builder->AddFESpace(std::string("Vector_L2"), std::string("L2_3D_P0"), 3); + problem_builder->AddGridFunction(std::string("magnetic_vector_potential"), std::string("HCurl")); + problem_builder->AddGridFunction(std::string("source_grad_phi"), std::string("HCurl")); + problem_builder->AddGridFunction(std::string("magnetic_flux_density"), std::string("HDiv")); + + problem_builder->AddGridFunction(std::string("magnetic_force"), std::string("Vector_L2")); + problem_builder->RegisterMagneticFluxDensityAux("magnetic_flux_density"); + + int outer_sphere_id = 101; + int inner_sphere_id = 102; + int sphere_volume_id = 100; + mfem::Array boundary_marker; + boundary_marker.Append(outer_sphere_id); + boundary_marker.Append(inner_sphere_id); + mfem::Array volume_marker; + volume_marker.Append(sphere_volume_id); + + hephaestus::Coefficients coefficients = defineCoefficients(); + problem_builder->SetCoefficients(coefficients); + + hephaestus::Sources sources = defineSources(); + problem_builder->SetSources(sources); + + hephaestus::Outputs outputs = defineOutputs(); + problem_builder->SetOutputs(outputs); + + auto magnetic_force_monitor = std::make_shared( + "magnetic_flux_density", "magnetic_vector_potential", boundary_marker, "magnetic_force" + // "magnetic_flux_density", "magnetic_vector_potential", volume_marker, "magnetic_force", + // false + ); + magnetic_force_monitor->SetPriority(2); + problem_builder->AddPostprocessor("MaxwellStressMonitor", magnetic_force_monitor); + + hephaestus::InputParameters solver_options; + solver_options.SetParam("Tolerance", float(1.0e-13)); + solver_options.SetParam("AbsTolerance", float(1.0e-16)); + solver_options.SetParam("MaxIter", (unsigned int)100); + problem_builder->SetSolverOptions(solver_options); + + problem_builder->FinalizeProblem(); + + auto problem = problem_builder->ReturnProblem(); + hephaestus::InputParameters exec_params; + exec_params.SetParam("VisualisationSteps", int(1)); + exec_params.SetParam("UseGLVis", true); + exec_params.SetParam("Problem", static_cast(problem.get())); + + auto executioner = std::make_unique(exec_params); + + hephaestus::logger.info("Created executioner"); + executioner->Execute(); + + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + double t; + double force; + int block_id; + for (std::size_t i = 0; i < magnetic_force_monitor->_times.Size(); ++i) + { + if (rank == 0) + { + t = magnetic_force_monitor->_times[i]; + block_id = magnetic_force_monitor->_attr_out[i]; + force = magnetic_force_monitor->_forces[i]; + hephaestus::logger.info("t = {} s, block {}, F = {} N", t, block_id, force); + } + } + + MPI_Finalize(); +} \ No newline at end of file diff --git a/src/auxsolvers/auxsolvers.hpp b/src/auxsolvers/auxsolvers.hpp index 94541380f..b44825632 100644 --- a/src/auxsolvers/auxsolvers.hpp +++ b/src/auxsolvers/auxsolvers.hpp @@ -11,6 +11,7 @@ #include "vector_gridfunction_dot_product_aux.hpp" #include "line_sampler_aux.hpp" #include "flux_monitor_aux.hpp" +#include "magnetic_force_aux.hpp" // Specify classes that perform auxiliary calculations on GridFunctions or // Coefficients. diff --git a/src/auxsolvers/magnetic_force_aux.cpp b/src/auxsolvers/magnetic_force_aux.cpp new file mode 100644 index 000000000..56af068b2 --- /dev/null +++ b/src/auxsolvers/magnetic_force_aux.cpp @@ -0,0 +1,188 @@ +#include "magnetic_force_aux.hpp" + +namespace hephaestus +{ +double +calcSurfaceForceDensity(mfem::ParGridFunction * b_field, + mfem::ParGridFunction * h_field, + int attr, + mfem::ParGridFunction & gf, + mfem::Coefficient & mu, + bool use_face_attr) +{ + double force = 0.0; + double total_force = 0.0; + double force_density = 0.0; + double area = 0.0; + + double air_permeability = M_PI * 4.0e-7; + + mfem::ParFiniteElementSpace * gf_fes = gf.ParFESpace(); + mfem::ParFiniteElementSpace * b_fes = b_field->ParFESpace(); + mfem::ParFiniteElementSpace * h_fes = h_field->ParFESpace(); + + mfem::ParMesh * mesh = gf_fes->GetParMesh(); + + mfem::Vector normal_vec, unit_normal_vec; + mfem::Array g_dof_ids; + + mfem::FaceElementTransformations * f_tr = nullptr; + + for (int i = 0; i < mesh->GetNBE(); i++) + { + f_tr = mesh->GetFaceElementTransformations(mesh->GetBdrElementFaceIndex(i)); + + if (use_face_attr) + { + if (mesh->GetBdrAttribute(i) != attr) + continue; + } + else + { + if (mesh->GetAttribute(f_tr->Elem1No) != attr) + { + continue; + } + } + const mfem::FiniteElement & elem = *b_fes->GetBE(i); + const mfem::IntegrationRule * ir = nullptr; + if (ir == nullptr) + { + const int order = 2 * elem.GetOrder() + 3; + ir = &mfem::IntRules.Get(f_tr->FaceGeom, order); + } + const int space_dim = 3; + + normal_vec.SetSize(space_dim); + unit_normal_vec.SetSize(space_dim); + + const mfem::FiniteElement * el = gf_fes->GetFE(i); + + double force_i = 0.0; + double area_i = 0.0; + double force_density_i = 0.0; + double force_density_j = 0.0; + + mfem::Vector force_vec_elem(space_dim); + force_vec_elem *= 0.0; + + // get dofs for writing to gridfunction + // gf_fes->GetBdrElementVDofs(i, g_dof_ids); + gf_fes->GetElementVDofs(f_tr->Elem1No, g_dof_ids); + + for (int j = 0; j < ir->GetNPoints(); j++) + { + mfem::Vector force_vec_ip(space_dim); + const mfem::IntegrationPoint & ip = ir->IntPoint(j); + double face_weight(0.0); + mfem::IntegrationPoint eip; + f_tr->Loc1.Transform(ip, eip); + f_tr->Face->SetIntPoint(&ip); + mfem::CalcOrtho(f_tr->Face->Jacobian(), normal_vec); + face_weight = f_tr->Face->Weight(); + f_tr->Elem1->SetIntPoint(&eip); + + // setup empty vectors + mfem::Vector b_vec(space_dim); + mfem::Vector h_vec(space_dim); + mfem::Vector h_tang(space_dim); + + // get vector values at integration point + b_field->GetVectorValue(*f_tr->Elem1, ip, b_vec); + h_field->GetVectorValue(*f_tr->Elem1, ip, h_vec); + double sphere_permeability = mu.Eval(*f_tr->Elem1, ip); + + // compute b normal component + unit_normal_vec.Set(1.0 / face_weight, normal_vec); + double b_normal_val = b_vec * unit_normal_vec; + double h_normal_val = h_vec * unit_normal_vec; + for (int k = 0; k < space_dim; ++k) + { + h_tang(k) = h_vec(k) - (unit_normal_vec(k) * h_normal_val); + } + + double term_1(0.0); + double term_2(0.0); + term_1 = (b_normal_val * b_normal_val) * (1.0 / air_permeability - 1.0 / sphere_permeability); + term_2 = (h_tang * h_tang) * (air_permeability - sphere_permeability); + + // Measure the area of the boundary + area += ip.weight * face_weight; + + force_density_j = ((term_1 - term_2) / 2.0); + // take y component of force for hollow sphere/levitation example + force_vec_ip.Set(((term_1 - term_2) / 2.0) * ip.weight * face_weight, unit_normal_vec); + double force_j = force_vec_ip(1); + force_vec_elem += force_vec_ip; + + force_i += force_j; + force_density_i += force_density_j; + } + force_density += force_density_i; + // write force, works for Vector_L2. Does not work with Vector_H1. + for (int j = 0; j < g_dof_ids.Size(); j++) + { + int ldof = g_dof_ids[j]; + gf(ldof) = force_vec_elem(j); + } + force += force_i; + } + + MPI_Allreduce(&force, &total_force, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + + return total_force; +} + +// ************************************************************************** // + +MagneticForceAux::MagneticForceAux(std::string b_name, + std::string h_name, + mfem::Array attr, + std::string coef_name, + bool use_face_attr) + : _b_name(std::move(b_name)), + _h_name(std::move(h_name)), + _coef_name(std::move(coef_name)), + _attr(std::move(attr)), + _use_face_attr(use_face_attr) +{ +} + +void +MagneticForceAux::Init(const hephaestus::GridFunctions & gridfunctions, + hephaestus::Coefficients & coefficients) +{ + _b_gf = gridfunctions.Get(_b_name); + _h_gf = gridfunctions.Get(_h_name); + _mu_coef = coefficients._scalars.Get("magnetic_permeability"); + + if (gridfunctions.Has(_coef_name)) + { + _gf = gridfunctions.Get(_coef_name); + } + // init field + *_gf = 0.0; +} + +// ************************************************************************** // + +void +MagneticForceAux::Solve(double t) +{ + double force; + + if (_gf != nullptr) + { + for (int i = 0; i < _attr.Size(); ++i) + { + force = calcSurfaceForceDensity(_b_gf, _h_gf, _attr[i], *_gf, *_mu_coef, _use_face_attr); + _times.Append(t); + _attr_out.Append(_attr[i]); + _forces.Append(force); + } + } +} + +// ************************************************************************** // + +} // namespace hephaestus diff --git a/src/auxsolvers/magnetic_force_aux.hpp b/src/auxsolvers/magnetic_force_aux.hpp new file mode 100644 index 000000000..00a03f7a0 --- /dev/null +++ b/src/auxsolvers/magnetic_force_aux.hpp @@ -0,0 +1,54 @@ +#pragma once +#include "auxsolver_base.hpp" + +// Specify postprocessors that depend on one or more gridfunctions +namespace hephaestus +{ + +double calcSurfaceForceDensity(mfem::ParGridFunction * b_field, + mfem::ParGridFunction * h_field, + int attr, + mfem::Coefficient & q, + mfem::Coefficient & mu, + bool use_face_attr); + +// Class to calculate and store the flux of a vector GridFunction through a surface +// at each timestep, optionally scaled by a coefficient. +class MagneticForceAux : public AuxSolver +{ + +public: + MagneticForceAux() = default; + MagneticForceAux(std::string b_name, + std::string h_name, + mfem::Array attr, + std::string coef_name = "", + bool use_face_attr = true); + + ~MagneticForceAux() override = default; + + void Init(const hephaestus::GridFunctions & gridfunctions, + hephaestus::Coefficients & coefficients) override; + + void Solve(double t = 0.0) override; + + // void WriteForces(std::string fname, mfem::ParGridFunction & gf, int attr); + + std::string _b_name; // name of the vector variable + std::string _h_name; // name of the vector variable + std::string _coef_name; // name of the coefficient + + mfem::Array _times; + mfem::Array _attr_out; + mfem::Array _forces; + + mfem::ParGridFunction * _b_gf{nullptr}; + mfem::ParGridFunction * _h_gf{nullptr}; + mfem::ParGridFunction * _gf{nullptr}; + mfem::Coefficient * _mu_coef{nullptr}; + + mfem::Array _attr; + bool _use_face_attr; +}; + +} // namespace hephaestus