diff --git a/libmeanfield/impl/operators/contexts/gravity_field_context.cpp b/libmeanfield/impl/operators/contexts/gravity_field_context.cpp index da18b40..0ec2e9d 100644 --- a/libmeanfield/impl/operators/contexts/gravity_field_context.cpp +++ b/libmeanfield/impl/operators/contexts/gravity_field_context.cpp @@ -7,72 +7,55 @@ module mean_field; import :operators.context.gravity_field; namespace { + using DomainSchema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema; + void validate_displacement( - const mean_field::fem::FEM &f, - const mfem::Vector &displacement_true + const mean_field::field::FieldDofMap &displacement_map, + const mfem::Vector &displacement ) { MFEM_VERIFY( - f.displacementFes != nullptr, "GravityFieldGeometryContext requires the " - "displacement finite-element space." - ); - MFEM_VERIFY( - displacement_true.Size() == f.displacementFes->GetTrueVSize(), + displacement.Size() == displacement_map.reduced_size(), "GravityFieldGeometryContext received a displacement vector with " "the " "wrong size." ); - for (int i = 0; i < displacement_true.Size(); ++i) { + for (int i = 0; i < displacement.Size(); ++i) { MFEM_VERIFY( - std::isfinite(displacement_true(i)), "GravityFieldGeometryContext received a non-finite " - "displacement " - "value." + std::isfinite(displacement(i)), "GravityFieldGeometryContext received a non-finite " + "displacement " + "value." ); } } void validate_linearization_state( - const mean_field::fem::FEM &f, + const mean_field::field::FieldDofMap &density_map, + const mean_field::field::FieldDofMap &displacement_map, + const mean_field::field::FieldDofMap &gravity_gradient_map, + const mean_field::field::FieldDofMap &gravity_potential_map, const mean_field::operators::context::gravity_field::GravityFieldStateView &state ) { MFEM_VERIFY( - f.densityFes != nullptr, "GravityFieldLinearizationContext " - "requires the density finite-element " - "space." - ); - MFEM_VERIFY( - f.gravityPotentialFes != nullptr, "GravityFieldLinearizationContext requires the gravity-potential " - "finite-element space." - ); - MFEM_VERIFY( - f.gravityFluxFes != nullptr, "GravityFieldLinearizationContext requires the " - "gravity-gradient finite-element space." - ); - MFEM_VERIFY( - f.displacementFes != nullptr, "GravityFieldLinearizationContext requires " - "the displacement finite-element space." - ); - - MFEM_VERIFY( - state.density.Size() == f.densityFes->GetTrueVSize(), + state.density.Size() == density_map.reduced_size(), "GravityFieldLinearizationContext received a density vector with " "the " "wrong size." ); MFEM_VERIFY( - state.displacement.Size() == f.displacementFes->GetTrueVSize(), + state.displacement.Size() == displacement_map.reduced_size(), "GravityFieldLinearizationContext received a displacement vector " "with " "the wrong size." ); MFEM_VERIFY( - state.gravity_gradient.Size() == f.gravityFluxFes->GetTrueVSize(), + state.gravity_gradient.Size() == gravity_gradient_map.reduced_size(), "GravityFieldLinearizationContext received a gravity-gradient " "vector " "with the wrong size." ); MFEM_VERIFY( - state.gravity_potential.Size() == f.gravityPotentialFes->GetTrueVSize(), + state.gravity_potential.Size() == gravity_potential_map.reduced_size(), "GravityFieldLinearizationContext received a gravity-potential " "vector " "with the wrong size." @@ -116,7 +99,12 @@ namespace mean_field::operators::context::gravity_field { const mapping::DomainMapperStateless &domain_mapper ) : m_fem(f), - m_domain_mapper(domain_mapper) { + m_domain_mapper(domain_mapper), + m_displacement_map( + field::make_field_dof_map< + field::Displacement, + DomainSchema>(*f.displacementFes) + ) { MFEM_VERIFY(f.mesh != nullptr, "GravityFieldGeometryContext requires a mesh."); MFEM_VERIFY( f.gravityFluxFes != nullptr, "GravityFieldGeometryContext requires the " @@ -153,11 +141,11 @@ namespace mean_field::operators::context::gravity_field { } GravityFieldGeometryPreparation GravityFieldGeometryContext::Prepare( - const mfem::Vector &displacement_true, + const mfem::Vector &displacement, const DiscretizationRevision discretization_revision, const DisplacementRevision displacement_revision ) { - validate_displacement(m_fem, displacement_true); + validate_displacement(m_displacement_map, displacement); if (m_is_prepared) { MFEM_VERIFY( @@ -185,8 +173,8 @@ namespace mean_field::operators::context::gravity_field { auto mass_operator = std::make_unique(m_fem, m_domain_mapper); auto source_operator = std::make_unique(m_fem, m_domain_mapper); - mass_operator->Prepare(displacement_true); - source_operator->Prepare(displacement_true); + mass_operator->Prepare(displacement); + source_operator->Prepare(displacement); m_mass_operator = std::move(mass_operator); m_source_operator = std::move(source_operator); @@ -204,14 +192,15 @@ namespace mean_field::operators::context::gravity_field { "operator." ); - m_mass_operator->Prepare(displacement_true); - m_source_operator->Prepare(displacement_true); + m_mass_operator->Prepare(displacement); + m_source_operator->Prepare(displacement); preparation.rebuilt_mass_operator = true; preparation.rebuilt_source_operator = true; } - m_displacement_true = displacement_true; + m_displacement_true.SetSize(m_displacement_map.full_size()); + m_displacement_map.scatter(displacement, m_displacement_true); m_discretization_revision = discretization_revision; m_displacement_revision = displacement_revision; m_is_prepared = true; @@ -242,7 +231,7 @@ namespace mean_field::operators::context::gravity_field { return *m_source_operator; } - const mfem::Vector &GravityFieldGeometryContext::GetDisplacement() const { + const mfem::Vector &GravityFieldGeometryContext::GetDisplacementTrue() const { MFEM_VERIFY( m_is_prepared, "GravityFieldGeometryContext must be prepared before " "accessing its displacement." @@ -250,6 +239,10 @@ namespace mean_field::operators::context::gravity_field { return m_displacement_true; } + const field::FieldDofMap &GravityFieldGeometryContext::GetDisplacementMap() const noexcept { + return m_displacement_map; + } + DiscretizationRevision GravityFieldGeometryContext::GetDiscretizationRevision() const noexcept { return m_discretization_revision; } @@ -270,6 +263,21 @@ namespace mean_field::operators::context::gravity_field { m_geometry_context( f, domain_mapper + ), + m_density_map( + field::make_field_dof_map< + field::Density, + DomainSchema>(*f.densityFes) + ), + m_gravity_gradient_map( + field::make_field_dof_map< + field::Gravity, + DomainSchema>(*f.gravityFluxFes) + ), + m_gravity_potential_map( + field::make_field_dof_map< + field::Gravity, + DomainSchema>(*f.gravityPotentialFes) ) { MFEM_VERIFY( f.densityFes != nullptr, "GravityFieldLinearizationContext " @@ -294,7 +302,10 @@ namespace mean_field::operators::context::gravity_field { const GravityFieldStateView &state, const GravityFieldRevisions &revisions ) { - validate_linearization_state(m_fem, state); + validate_linearization_state( + m_density_map, m_geometry_context.GetDisplacementMap(), m_gravity_gradient_map, m_gravity_potential_map, + state + ); if (m_is_prepared) { MFEM_VERIFY( @@ -338,12 +349,14 @@ namespace mean_field::operators::context::gravity_field { m_geometry_context.Prepare(state.displacement, revisions.discretization, revisions.displacement); if (density_changed) { - m_density_true = state.density; + m_density_true.SetSize(m_density_map.full_size()); + m_density_map.scatter(state.density, m_density_true); report.updated_density = true; } if (gravity_gradient_changed) { - m_gravity_gradient_true = state.gravity_gradient; + m_gravity_gradient_true.SetSize(m_gravity_gradient_map.full_size()); + m_gravity_gradient_map.scatter(state.gravity_gradient, m_gravity_gradient_true); report.updated_gravity_gradient = true; } @@ -361,7 +374,7 @@ namespace mean_field::operators::context::gravity_field { return m_geometry_context; } - const mfem::Vector &GravityFieldLinearizationContext::GetDensity() const { + const mfem::Vector &GravityFieldLinearizationContext::GetDensityTrue() const { MFEM_VERIFY( m_is_prepared, "GravityFieldLinearizationContext must be prepared " "before accessing its density." @@ -369,7 +382,7 @@ namespace mean_field::operators::context::gravity_field { return m_density_true; } - const mfem::Vector &GravityFieldLinearizationContext::GetGravityGradient() const { + const mfem::Vector &GravityFieldLinearizationContext::GetGravityGradientTrue() const { MFEM_VERIFY( m_is_prepared, "GravityFieldLinearizationContext must be prepared " "before accessing its gravity gradient." @@ -377,6 +390,22 @@ namespace mean_field::operators::context::gravity_field { return m_gravity_gradient_true; } + const field::FieldDofMap &GravityFieldLinearizationContext::GetDensityMap() const noexcept { + return m_density_map; + } + + const field::FieldDofMap &GravityFieldLinearizationContext::GetDisplacementMap() const noexcept { + return m_geometry_context.GetDisplacementMap(); + } + + const field::FieldDofMap &GravityFieldLinearizationContext::GetGravityGradientMap() const noexcept { + return m_gravity_gradient_map; + } + + const field::FieldDofMap &GravityFieldLinearizationContext::GetGravityPotentialMap() const noexcept { + return m_gravity_potential_map; + } + const GravityFieldRevisions &GravityFieldLinearizationContext::GetRevisions() const { MFEM_VERIFY( m_is_prepared, "GravityFieldLinearizationContext must be prepared " diff --git a/libmeanfield/impl/operators/contexts/hydrostatic_equilibrium_context.cpp b/libmeanfield/impl/operators/contexts/hydrostatic_equilibrium_context.cpp index 55ac617..0ebca61 100644 --- a/libmeanfield/impl/operators/contexts/hydrostatic_equilibrium_context.cpp +++ b/libmeanfield/impl/operators/contexts/hydrostatic_equilibrium_context.cpp @@ -9,6 +9,8 @@ module mean_field; import :operators.context.hydrostatic_equilibrium; namespace { + using DomainSchema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema; + void validate_finite_vector( const mfem::Vector &vector, const char *message @@ -19,38 +21,24 @@ namespace { } void validate_state( - const mean_field::fem::FEM &f, + const mean_field::field::FieldDofMap &enthalpyMap, + const mean_field::field::FieldDofMap &gravityPotentialMap, + const mean_field::field::FieldDofMap &displacementMap, const mean_field::operators::context::hydrostatic::HydrostaticEquilibriumStateView &state ) { MFEM_VERIFY( - f.enthalpyFes != nullptr, "HydrostaticEquilibriumContext requires the " - "enthalpy finite-element space." + state.enthalpy.Size() == enthalpyMap.reduced_size(), "HydrostaticEquilibriumContext received a supported " + "enthalpy vector with the wrong size." ); MFEM_VERIFY( - f.gravityPotentialFes != nullptr, "HydrostaticEquilibriumContext requires the " - "gravity-potential finite-element space." + state.gravityPotential.Size() == gravityPotentialMap.reduced_size(), + "HydrostaticEquilibriumContext received a supported gravity-potential vector with the wrong size." ); MFEM_VERIFY( - f.displacementFes != nullptr, "HydrostaticEquilibriumContext requires the " - "displacement finite-element space." - ); - - MFEM_VERIFY( - state.enthalpy.Size() == f.enthalpyFes->GetTrueVSize(), "HydrostaticEquilibriumContext received an " - "enthalpy vector with the wrong size." - ); - - MFEM_VERIFY( - state.gravityPotential.Size() == f.gravityPotentialFes->GetTrueVSize(), - "HydrostaticEquilibriumContext received a " - "gravity-potential vector with the wrong size." - ); - - MFEM_VERIFY( - state.displacement.Size() == f.displacementFes->GetTrueVSize(), "HydrostaticEquilibriumContext received a " - "displacement vector with the wrong size." + state.displacement.Size() == displacementMap.reduced_size(), + "HydrostaticEquilibriumContext received a supported displacement vector with the wrong size." ); validate_finite_vector( @@ -90,7 +78,16 @@ namespace mean_field::operators::context::hydrostatic { const mapping::DomainMapperStateless &domainMapper ) : m_f(f), - m_domainMapper(domainMapper) { + m_domainMapper(domainMapper), + m_enthalpyMap( + field::make_field_dof_map(*f.enthalpyFes) + ), + m_gravityPotentialMap( + field::make_field_dof_map(*f.gravityPotentialFes) + ), + m_displacementMap( + field::make_field_dof_map(*f.displacementFes) + ) { MFEM_VERIFY(m_f.mesh != nullptr, "HydrostaticEquilibriumContext requires a mesh."); MFEM_VERIFY( @@ -113,13 +110,17 @@ namespace mean_field::operators::context::hydrostatic { "domain-mapper dimension does not match the mesh " "dimension." ); + + m_baseEnthalpyTrue.SetSize(m_enthalpyMap.full_size()); + m_baseGravityPotentialTrue.SetSize(m_gravityPotentialMap.full_size()); + m_displacementTrue.SetSize(m_displacementMap.full_size()); } HydrostaticPreparationReport HydrostaticEquilibriumContext::Prepare( const HydrostaticEquilibriumStateView &state, const HydrostaticEquilibriumDependencies &dependencies ) { - validate_state(m_f, state); + validate_state(m_enthalpyMap, m_gravityPotentialMap, m_displacementMap, state); if (m_isPrepared) { validate_dependency_transition( @@ -188,18 +189,18 @@ namespace mean_field::operators::context::hydrostatic { report.preparedBaseState = baseStatePreparationRequired; if (staticChanged || enthalpyChanged) { - m_baseEnthalpyTrue = state.enthalpy; + m_enthalpyMap.scatter(state.enthalpy, m_baseEnthalpyTrue); report.updatedEnthalpy = true; } if (staticChanged || gravityPotentialChanged) { - m_baseGravityPotentialTrue = state.gravityPotential; + m_gravityPotentialMap.scatter(state.gravityPotential, m_baseGravityPotentialTrue); report.updatedGravityPotential = true; } if (geometryPreparationRequired) { - m_displacementTrue = state.displacement; + m_displacementMap.scatter(state.displacement, m_displacementTrue); report.updatedDisplacement = true; } @@ -249,6 +250,18 @@ namespace mean_field::operators::context::hydrostatic { return m_statistics; } + const field::FieldDofMap &HydrostaticEquilibriumContext::GetEnthalpyMap() const noexcept { + return m_enthalpyMap; + } + + const field::FieldDofMap &HydrostaticEquilibriumContext::GetGravityPotentialMap() const noexcept { + return m_gravityPotentialMap; + } + + const field::FieldDofMap &HydrostaticEquilibriumContext::GetDisplacementMap() const noexcept { + return m_displacementMap; + } + const mfem::Vector &HydrostaticEquilibriumContext::GetBaseEnthalpyTrue() const { VerifyPrepared(); return m_baseEnthalpyTrue; diff --git a/libmeanfield/impl/operators/contexts/rotation_displacement_force_context.cpp b/libmeanfield/impl/operators/contexts/rotation_displacement_force_context.cpp index 43187e2..47a3e86 100644 --- a/libmeanfield/impl/operators/contexts/rotation_displacement_force_context.cpp +++ b/libmeanfield/impl/operators/contexts/rotation_displacement_force_context.cpp @@ -9,6 +9,8 @@ module mean_field; import :operators.context.rotational_displacement_force; namespace { + using DomainSchema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema; + void validate_finite_vector( const mfem::Vector &vector, const char *message @@ -33,7 +35,17 @@ namespace mean_field::operators::context::rotational_displacement_force { const fem::FEM &f, const mapping::DomainMapperStateless &domainMapper ) - : m_f(f) { + : m_f(f), + m_densityMap( + field::make_field_dof_map< + field::Density, + DomainSchema>(*f.densityFes) + ), + m_displacementMap( + field::make_field_dof_map< + field::Displacement, + DomainSchema>(*f.displacementFes) + ) { MFEM_VERIFY( m_f.mesh != nullptr, "RotationalDisplacementForceLinearizationContext requires a " "mesh." @@ -61,13 +73,13 @@ namespace mean_field::operators::context::rotational_displacement_force { const RotationalDisplacementForceDependencies &dependencies ) { MFEM_VERIFY( - state.density.Size() == m_f.densityFes->GetTrueVSize(), + state.density.Size() == m_densityMap.reduced_size(), "RotationalDisplacementForceLinearizationContext received a " "density vector with the wrong size." ); MFEM_VERIFY( - state.displacement.Size() == m_f.displacementFes->GetTrueVSize(), + state.displacement.Size() == m_displacementMap.reduced_size(), "RotationalDisplacementForceLinearizationContext received a " "displacement vector with the wrong size." ); @@ -132,12 +144,14 @@ namespace mean_field::operators::context::rotational_displacement_force { report.preparedBaseState = baseStatePreparationRequired; if (discretizationChanged || densityChanged) { - m_baseDensityTrue = state.density; + m_baseDensityTrue.SetSize(m_densityMap.full_size()); + m_densityMap.scatter(state.density, m_baseDensityTrue); report.updatedDensity = true; } if (geometryPreparationRequired) { - m_displacementTrue = state.displacement; + m_displacementTrue.SetSize(m_displacementMap.full_size()); + m_displacementMap.scatter(state.displacement, m_displacementTrue); report.updatedDisplacement = true; } @@ -194,6 +208,14 @@ namespace mean_field::operators::context::rotational_displacement_force { return m_displacementTrue; } + const field::FieldDofMap &RotationalDisplacementForceLinearizationContext::GetDensityMap() const noexcept { + return m_densityMap; + } + + const field::FieldDofMap &RotationalDisplacementForceLinearizationContext::GetDisplacementMap() const noexcept { + return m_displacementMap; + } + void RotationalDisplacementForceLinearizationContext::VerifyPrepared() const { MFEM_VERIFY( m_isPrepared, "RotationalDisplacementForceLinearizationContext has not been " diff --git a/libmeanfield/impl/operators/gravity_field.cpp b/libmeanfield/impl/operators/gravity_field.cpp index 62dfa34..e2b7141 100644 --- a/libmeanfield/impl/operators/gravity_field.cpp +++ b/libmeanfield/impl/operators/gravity_field.cpp @@ -12,18 +12,16 @@ import :operators.kernels.gravity_field; namespace { using namespace mean_field; - int get_state_width(const mfem::Array &state_true_offsets) { - MFEM_VERIFY(state_true_offsets.Size() >= 2, "The coupled state requires at least one block."); - MFEM_VERIFY(state_true_offsets[0] == 0, "The coupled state offsets must begin at zero."); + int get_state_width(const mfem::Array &state_offsets) { + MFEM_VERIFY(state_offsets.Size() >= 2, "The coupled state requires at least one block."); + MFEM_VERIFY(state_offsets[0] == 0, "The coupled state offsets must begin at zero."); - for (int i = 0; i < state_true_offsets.Size() - 1; ++i) { - MFEM_VERIFY( - state_true_offsets[i + 1] >= state_true_offsets[i], "The coupled state offsets must be nondecreasing." - ); + for (int i = 0; i < state_offsets.Size() - 1; ++i) { + MFEM_VERIFY(state_offsets[i + 1] >= state_offsets[i], "The coupled state offsets must be nondecreasing."); } - MFEM_VERIFY(state_true_offsets.Last() > 0, "The coupled state cannot be empty."); - return state_true_offsets.Last(); + MFEM_VERIFY(state_offsets.Last() > 0, "The coupled state cannot be empty."); + return state_offsets.Last(); } int get_gravity_residual_height(const fem::FEM &f) { @@ -37,29 +35,33 @@ namespace { "space (L2: Lebesgue " "space of square-integrable functions)." ); - return f.gravityFluxFes->GetTrueVSize() + f.gravityPotentialFes->GetTrueVSize(); + using DomainSchema = utils::domain::CoreEnvelopeVacuumDomainSchema; + return field::make_field_dof_map(*f.gravityFluxFes).reduced_size() + + field::make_field_dof_map(*f.gravityPotentialFes).reduced_size(); } mfem::Array make_gravity_residual_offsets(const fem::FEM &f) { - mfem::Array offsets(operators::gravity_residual_block_count + 1); - offsets[0] = 0; - offsets[1] = f.gravityFluxFes->GetTrueVSize(); - offsets[2] = offsets[1] + f.gravityPotentialFes->GetTrueVSize(); + mfem::Array offsets(utils::blocks::gravity_field_form::residual_block_count + 1); + offsets[0] = 0; + using DomainSchema = utils::domain::CoreEnvelopeVacuumDomainSchema; + offsets[1] = field::make_field_dof_map(*f.gravityFluxFes).reduced_size(); + offsets[2] = + offsets[1] + field::make_field_dof_map(*f.gravityPotentialFes).reduced_size(); return offsets; } template int get_state_block_size( - const mfem::Array &state_true_offsets, + const mfem::Array &state_offsets, const utils::blocks::value_block ) { - MFEM_VERIFY(index + 1 < state_true_offsets.Size(), "Value block is not present in the state offsets."); - return state_true_offsets[index + 1] - state_true_offsets[index]; + MFEM_VERIFY(index + 1 < state_offsets.Size(), "Value block is not present in the state offsets."); + return state_offsets[index + 1] - state_offsets[index]; } void validate_state_offsets( const fem::FEM &f, - const mfem::Array &state_true_offsets + const mfem::Array &state_offsets ) { MFEM_VERIFY(f.densityFes != nullptr, "GravityFieldOperator requires the density finite-element space."); MFEM_VERIFY( @@ -78,26 +80,32 @@ namespace { utils::blocks::get_value_block
(utils::blocks::gravity_field.poisson_term); MFEM_VERIFY( - state_true_offsets.Size() == form::value_block_count + 1, + state_offsets.Size() == form::value_block_count + 1, "The gravity state offsets do not match gravity_field_form." ); + using DomainSchema = utils::domain::CoreEnvelopeVacuumDomainSchema; + const auto density_map = field::make_field_dof_map(*f.densityFes); + const auto displacement_map = field::make_field_dof_map(*f.displacementFes); + const auto flux_map = field::make_field_dof_map(*f.gravityFluxFes); + const auto potential_map = field::make_field_dof_map(*f.gravityPotentialFes); + MFEM_VERIFY( - get_state_block_size(state_true_offsets, density_block) == f.densityFes->GetTrueVSize(), + get_state_block_size(state_offsets, density_block) == density_map.reduced_size(), "The density block does not match the density finite-element space." ); MFEM_VERIFY( - get_state_block_size(state_true_offsets, displacement_block) == f.displacementFes->GetTrueVSize(), + get_state_block_size(state_offsets, displacement_block) == displacement_map.reduced_size(), "The displacement block does not match the displacement " "finite-element " "space." ); MFEM_VERIFY( - get_state_block_size(state_true_offsets, gravity_gradient_block) == f.gravityFluxFes->GetTrueVSize(), + get_state_block_size(state_offsets, gravity_gradient_block) == flux_map.reduced_size(), "The gravity-gradient block does not match the RT finite-element " "space." ); MFEM_VERIFY( - get_state_block_size(state_true_offsets, gravity_potential_block) == f.gravityPotentialFes->GetTrueVSize(), + get_state_block_size(state_offsets, gravity_potential_block) == potential_map.reduced_size(), "The gravity-potential block does not match the potential " "finite-element space." ); @@ -162,18 +170,18 @@ namespace mean_field::operators { fem::FEM &f, const mapping::DomainMapperStateless &domain_mapper, context::gravity_field::GravityFieldLinearizationContext &linearization_context, - const mfem::Array &state_true_offsets, + const mfem::Array &state_offsets, GravityFieldJacobianOperator &jacobian ) : Operator( get_gravity_residual_height(f), - get_state_width(state_true_offsets) + get_state_width(state_offsets) ), m_fem(f), m_domain_mapper(domain_mapper), m_linearization_context(linearization_context), - m_state_true_offsets(state_true_offsets), - m_residual_true_offsets(make_gravity_residual_offsets(f)), + m_state_offsets(state_offsets), + m_residual_offsets(make_gravity_residual_offsets(f)), m_jacobian(jacobian) { MFEM_VERIFY(f.mesh != nullptr, "GravityFieldOperator requires a mesh."); MFEM_VERIFY( @@ -198,7 +206,7 @@ namespace mean_field::operators { "dimension." ); - validate_state_offsets(f, m_state_true_offsets); + validate_state_offsets(f, m_state_offsets); validate_gravity_context(f); bool has_vacuum_domain = false; @@ -212,11 +220,9 @@ namespace mean_field::operators { MFEM_VERIFY(has_vacuum_domain, "GravityFieldOperator requires a compactified vacuum domain."); MFEM_VERIFY( - m_residual_true_offsets.Last() == Height(), "The gravity residual offsets do not match the operator height." - ); - MFEM_VERIFY( - m_state_true_offsets.Last() == Width(), "The coupled state offsets do not match the operator width." + m_residual_offsets.Last() == Height(), "The gravity residual offsets do not match the operator height." ); + MFEM_VERIFY(m_state_offsets.Last() == Width(), "The coupled state offsets do not match the operator width."); } context::gravity_field::GravityFieldPreparationReport GravityFieldOperator::Prepare( @@ -238,12 +244,11 @@ namespace mean_field::operators { "preparation state with the wrong size." ); - const mfem::Vector density = make_read_only_value_view(state, m_state_true_offsets, density_block); - const mfem::Vector displacement = make_read_only_value_view(state, m_state_true_offsets, displacement_block); - const mfem::Vector gravity_gradient = - make_read_only_value_view(state, m_state_true_offsets, gravity_gradient_block); + const mfem::Vector density = make_read_only_value_view(state, m_state_offsets, density_block); + const mfem::Vector displacement = make_read_only_value_view(state, m_state_offsets, displacement_block); + const mfem::Vector gravity_gradient = make_read_only_value_view(state, m_state_offsets, gravity_gradient_block); const mfem::Vector gravity_potential = - make_read_only_value_view(state, m_state_true_offsets, gravity_potential_block); + make_read_only_value_view(state, m_state_offsets, gravity_potential_block); return m_linearization_context.Prepare( {.density = density, @@ -254,12 +259,12 @@ namespace mean_field::operators { ); } - const mfem::Array &GravityFieldOperator::GetStateTrueOffsets() const noexcept { - return m_state_true_offsets; + const mfem::Array &GravityFieldOperator::GetStateOffsets() const noexcept { + return m_state_offsets; } - const mfem::Array &GravityFieldOperator::GetResidualTrueOffsets() const noexcept { - return m_residual_true_offsets; + const mfem::Array &GravityFieldOperator::GetResidualOffsets() const noexcept { + return m_residual_offsets; } void GravityFieldOperator::ApplyGravityUnknowns( @@ -277,12 +282,12 @@ namespace mean_field::operators { MFEM_VERIFY(geometry_context.IsPrepared(), "GravityFieldOperator received an unprepared geometry context."); MFEM_VERIFY( - gravity_gradient.Size() == m_fem.gravityFluxFes->GetTrueVSize(), + gravity_gradient.Size() == geometry_context.GetMassOperator().GetFluxMap().reduced_size(), "GravityFieldOperator received a gravity-gradient vector with the " "wrong size." ); MFEM_VERIFY( - gravity_potential.Size() == m_fem.gravityPotentialFes->GetTrueVSize(), + gravity_potential.Size() == geometry_context.GetSourceOperator().GetPotentialMap().reduced_size(), "GravityFieldOperator received a gravity-potential vector with the " "wrong size." ); @@ -291,15 +296,25 @@ namespace mean_field::operators { action = 0.0; mfem::Vector gravity_gradient_action = - make_residual_view(action, m_residual_true_offsets, gravity_gradient_residual_block); + make_residual_view(action, m_residual_offsets, gravity_gradient_residual_block); mfem::Vector gravity_poisson_action = - make_residual_view(action, m_residual_true_offsets, gravity_poisson_residual_block); - mfem::Vector transpose_divergence_action(gravity_gradient_action.Size()); + make_residual_view(action, m_residual_offsets, gravity_poisson_residual_block); + const field::FieldDofMap &flux_map = geometry_context.GetMassOperator().GetFluxMap(); + const field::FieldDofMap &potential_map = geometry_context.GetSourceOperator().GetPotentialMap(); + mfem::Vector potential_true(potential_map.full_size()); + mfem::Vector transpose_divergence_action_true(flux_map.full_size()); + mfem::Vector transpose_divergence_action(flux_map.reduced_size()); + mfem::Vector gradient_true(flux_map.full_size()); + mfem::Vector divergence_action_true(potential_map.full_size()); geometry_context.GetMassOperator().Mult(gravity_gradient, gravity_gradient_action); - m_fem.gravityContext.BT->Mult(gravity_potential, transpose_divergence_action); + potential_map.scatter(gravity_potential, potential_true); + m_fem.gravityContext.BT->Mult(potential_true, transpose_divergence_action_true); + flux_map.gather(transpose_divergence_action_true, transpose_divergence_action); gravity_gradient_action += transpose_divergence_action; - m_fem.gravityContext.b_form->Mult(gravity_gradient, gravity_poisson_action); + flux_map.scatter(gravity_gradient, gradient_true); + m_fem.gravityContext.b_form->Mult(gradient_true, divergence_action_true); + potential_map.gather(divergence_action_true, gravity_poisson_action); } void GravityFieldOperator::ApplyDensitySource( @@ -314,7 +329,7 @@ namespace mean_field::operators { MFEM_VERIFY(geometry_context.IsPrepared(), "GravityFieldOperator received an unprepared geometry context."); MFEM_VERIFY( - density.Size() == m_fem.densityFes->GetTrueVSize(), + density.Size() == geometry_context.GetSourceOperator().GetDensityMap().reduced_size(), "GravityFieldOperator received a density vector with the wrong " "size." ); @@ -323,7 +338,7 @@ namespace mean_field::operators { action = 0.0; mfem::Vector gravity_poisson_action = - make_residual_view(action, m_residual_true_offsets, gravity_poisson_residual_block); + make_residual_view(action, m_residual_offsets, gravity_poisson_residual_block); geometry_context.GetSourceOperator().Mult(density, gravity_poisson_action); } @@ -346,11 +361,10 @@ namespace mean_field::operators { m_linearization_context.IsPrepared(), "GravityFieldOperator must be prepared before Mult is called." ); - const mfem::Vector density = make_read_only_value_view(state, m_state_true_offsets, density_block); - const mfem::Vector gravity_gradient = - make_read_only_value_view(state, m_state_true_offsets, gravity_gradient_block); + const mfem::Vector density = make_read_only_value_view(state, m_state_offsets, density_block); + const mfem::Vector gravity_gradient = make_read_only_value_view(state, m_state_offsets, gravity_gradient_block); const mfem::Vector gravity_potential = - make_read_only_value_view(state, m_state_true_offsets, gravity_potential_block); + make_read_only_value_view(state, m_state_offsets, gravity_potential_block); const context::gravity_field::GravityFieldGeometryContext &geometry_context = m_linearization_context.GetGeometryContext(); @@ -393,7 +407,7 @@ namespace mean_field::operators { gravity_field_operator.Height() ), m_gravity_field_operator(gravity_field_operator), - m_gravity_true_offsets(gravity_field_operator.GetResidualTrueOffsets()), + m_gravity_offsets(gravity_field_operator.GetResidualOffsets()), m_gravity_field_geometry_context(gravity_field_geometry_context) { using form = utils::blocks::gravity_field_form; @@ -406,52 +420,48 @@ namespace mean_field::operators { constexpr auto gravity_poisson_residual_block = utils::blocks::get_residual_block(utils::blocks::gravity_field.poisson_term); - const mfem::Array &state_offsets = m_gravity_field_operator.GetStateTrueOffsets(); + const mfem::Array &state_offsets = m_gravity_field_operator.GetStateOffsets(); MFEM_VERIFY( state_offsets.Size() == form::value_block_count + 1, - "ReducedGravityFieldOperator received an invalid full-state layout." + "ReducedGravityFieldOperator received an invalid coupled-state layout." ); MFEM_VERIFY( - m_gravity_true_offsets.Size() == form::residual_block_count + 1, + m_gravity_offsets.Size() == form::residual_block_count + 1, "ReducedGravityFieldOperator received an invalid gravity-residual " "layout." ); - MFEM_VERIFY(state_offsets[0] == 0, "The full-state offsets must begin at zero."); - MFEM_VERIFY(m_gravity_true_offsets[0] == 0, "The reduced gravity offsets must begin at zero."); + MFEM_VERIFY(state_offsets[0] == 0, "The coupled-state offsets must begin at zero."); + MFEM_VERIFY(m_gravity_offsets[0] == 0, "The reduced gravity offsets must begin at zero."); MFEM_VERIFY( state_offsets.Last() == m_gravity_field_operator.Width(), - "The full-state offsets do not match the gravity-field operator " + "The coupled-state offsets do not match the gravity-field operator " "width." ); MFEM_VERIFY( - m_gravity_true_offsets.Last() == m_gravity_field_operator.Height(), + m_gravity_offsets.Last() == m_gravity_field_operator.Height(), "The reduced gravity offsets do not match the gravity-field " "operator " "height." ); MFEM_VERIFY(Width() == Height(), "ReducedGravityFieldOperator must be square."); - const int full_gradient_size = + const int state_gradient_size = state_offsets[static_cast(gravity_gradient_block) + 1] - state_offsets[gravity_gradient_block]; - const int full_potential_size = + const int state_potential_size = state_offsets[static_cast(gravity_potential_block) + 1] - state_offsets[gravity_potential_block]; - const int reduced_gradient_size = - m_gravity_true_offsets[static_cast(gravity_gradient_residual_block) + 1] - - m_gravity_true_offsets[gravity_gradient_residual_block]; - const int reduced_potential_size = - m_gravity_true_offsets[static_cast(gravity_poisson_residual_block) + 1] - - m_gravity_true_offsets[gravity_poisson_residual_block]; + const int gravity_gradient_size = m_gravity_offsets[static_cast(gravity_gradient_residual_block) + 1] - + m_gravity_offsets[gravity_gradient_residual_block]; + const int gravity_potential_size = m_gravity_offsets[static_cast(gravity_poisson_residual_block) + 1] - + m_gravity_offsets[gravity_poisson_residual_block]; MFEM_VERIFY( - full_gradient_size == reduced_gradient_size, - "The reduced gravity-gradient block does not match the full-state " - "gravity-gradient block." + state_gradient_size == gravity_gradient_size, "The gravity-gradient block does not match the coupled-state " + "gravity-gradient block." ); MFEM_VERIFY( - full_potential_size == reduced_potential_size, - "The reduced gravity-potential block does not match the Poisson " - "residual block." + state_potential_size == gravity_potential_size, "The gravity-potential block does not match the Poisson " + "residual block." ); SetDisplacement(displacement); @@ -475,10 +485,11 @@ namespace mean_field::operators { } m_gravity_field_geometry_context.Prepare(displacement, discretization_revision, displacement_revision); + m_displacement = displacement; } const mfem::Vector &ReducedGravityFieldOperator::GetDisplacement() const { - return m_gravity_field_geometry_context.GetDisplacement(); + return m_displacement; } void ReducedGravityFieldOperator::BuildRightHandSide( @@ -511,13 +522,13 @@ namespace mean_field::operators { ValidateGravityState(gravity_state); - const mfem::Vector gravity_gradient_true = - make_read_only_residual_view(gravity_state, m_gravity_true_offsets, gravity_gradient_residual_block); - const mfem::Vector gravity_potential_true = - make_read_only_residual_view(gravity_state, m_gravity_true_offsets, gravity_poisson_residual_block); + const mfem::Vector gravity_gradient = + make_read_only_residual_view(gravity_state, m_gravity_offsets, gravity_gradient_residual_block); + const mfem::Vector gravity_potential = + make_read_only_residual_view(gravity_state, m_gravity_offsets, gravity_poisson_residual_block); m_gravity_field_operator.ApplyGravityUnknowns( - gravity_gradient_true, gravity_potential_true, m_gravity_field_geometry_context, action + gravity_gradient, gravity_potential, m_gravity_field_geometry_context, action ); MFEM_VERIFY( @@ -543,8 +554,8 @@ namespace mean_field::operators { return m_gravity_field_geometry_context; } - const mfem::Array &ReducedGravityFieldOperator::GetGravityTrueOffsets() const noexcept { - return m_gravity_true_offsets; + const mfem::Array &ReducedGravityFieldOperator::GetGravityOffsets() const noexcept { + return m_gravity_offsets; } void ReducedGravityFieldOperator::ValidateDisplacement(const mfem::Vector &displacement) const { @@ -553,7 +564,7 @@ namespace mean_field::operators { constexpr auto displacement_block = utils::blocks::get_value_block(utils::blocks::displacement_field.geometry_term); - const mfem::Array &state_offsets = m_gravity_field_operator.GetStateTrueOffsets(); + const mfem::Array &state_offsets = m_gravity_field_operator.GetStateOffsets(); const int expected_size = state_offsets[static_cast(displacement_block) + 1] - state_offsets[displacement_block]; @@ -577,7 +588,7 @@ namespace mean_field::operators { constexpr auto density_block = utils::blocks::get_value_block(utils::blocks::density_field.mass_term); - const mfem::Array &state_offsets = m_gravity_field_operator.GetStateTrueOffsets(); + const mfem::Array &state_offsets = m_gravity_field_operator.GetStateOffsets(); const int expected_size = state_offsets[static_cast(density_block) + 1] - state_offsets[density_block]; MFEM_VERIFY( diff --git a/libmeanfield/impl/operators/gravity_field_jacobian.cpp b/libmeanfield/impl/operators/gravity_field_jacobian.cpp index 9242324..cc166f7 100644 --- a/libmeanfield/impl/operators/gravity_field_jacobian.cpp +++ b/libmeanfield/impl/operators/gravity_field_jacobian.cpp @@ -89,28 +89,38 @@ namespace { "form." ); + using DomainSchema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema; + const auto density_map = + mean_field::field::make_field_dof_map(*f.densityFes); + const auto displacement_map = + mean_field::field::make_field_dof_map(*f.displacementFes); + const auto flux_map = + mean_field::field::make_field_dof_map(*f.gravityFluxFes); + const auto potential_map = + mean_field::field::make_field_dof_map(*f.gravityPotentialFes); + MFEM_VERIFY( - get_block_size(state_offsets, density_block) == f.densityFes->GetTrueVSize(), + get_block_size(state_offsets, density_block) == density_map.reduced_size(), "The Jacobian density block has the wrong size." ); MFEM_VERIFY( - get_block_size(state_offsets, displacement_block) == f.displacementFes->GetTrueVSize(), + get_block_size(state_offsets, displacement_block) == displacement_map.reduced_size(), "The Jacobian displacement block has the wrong size." ); MFEM_VERIFY( - get_block_size(state_offsets, gravity_gradient_block) == f.gravityFluxFes->GetTrueVSize(), + get_block_size(state_offsets, gravity_gradient_block) == flux_map.reduced_size(), "The Jacobian gravity-gradient block has the wrong size." ); MFEM_VERIFY( - get_block_size(state_offsets, gravity_potential_block) == f.gravityPotentialFes->GetTrueVSize(), + get_block_size(state_offsets, gravity_potential_block) == potential_map.reduced_size(), "The Jacobian gravity-potential block has the wrong size." ); MFEM_VERIFY( - get_block_size(residual_offsets, gravity_gradient_residual_block) == f.gravityFluxFes->GetTrueVSize(), + get_block_size(residual_offsets, gravity_gradient_residual_block) == flux_map.reduced_size(), "The Jacobian gradient-residual block has the wrong size." ); MFEM_VERIFY( - get_block_size(residual_offsets, gravity_poisson_residual_block) == f.gravityPotentialFes->GetTrueVSize(), + get_block_size(residual_offsets, gravity_poisson_residual_block) == potential_map.reduced_size(), "The Jacobian Poisson-residual block has the wrong size." ); } @@ -121,18 +131,18 @@ namespace mean_field::operators { fem::FEM &f, const mapping::DomainMapperStateless &domain_mapper, const context::gravity_field::GravityFieldLinearizationContext &linearization_context, - const mfem::Array &state_true_offsets, - const mfem::Array &residual_true_offsets + const mfem::Array &state_offsets, + const mfem::Array &residual_offsets ) : Operator( - residual_true_offsets.Last(), - state_true_offsets.Last() + residual_offsets.Last(), + state_offsets.Last() ), m_fem(f), m_domain_mapper(domain_mapper), m_linearization_context(linearization_context), - m_state_true_offsets(state_true_offsets), - m_residual_true_offsets(residual_true_offsets) { + m_state_offsets(state_offsets), + m_residual_offsets(residual_offsets) { MFEM_VERIFY( f.densityFes != nullptr, "GravityFieldJacobianOperator requires the density finite-element " "space." @@ -166,7 +176,7 @@ namespace mean_field::operators { "dimension." ); - validate_layout(f, m_state_true_offsets, m_residual_true_offsets); + validate_layout(f, m_state_offsets, m_residual_offsets); } void GravityFieldJacobianOperator::Mult( @@ -198,50 +208,70 @@ namespace mean_field::operators { const context::gravity_field::GravityFieldGeometryContext &geometry_context = m_linearization_context.GetGeometryContext(); - const mfem::Vector &density = m_linearization_context.GetDensity(); - const mfem::Vector &displacement = geometry_context.GetDisplacement(); - const mfem::Vector &gravity_gradient = m_linearization_context.GetGravityGradient(); + const mfem::Vector &density = m_linearization_context.GetDensityTrue(); + const mfem::Vector &displacement = geometry_context.GetDisplacementTrue(); + const mfem::Vector &gravity_gradient = m_linearization_context.GetGravityGradientTrue(); - const mfem::Vector density_direction = - make_read_only_value_view(direction, m_state_true_offsets, density_block); + const mfem::Vector density_direction = make_read_only_value_view(direction, m_state_offsets, density_block); const mfem::Vector displacement_direction = - make_read_only_value_view(direction, m_state_true_offsets, displacement_block); + make_read_only_value_view(direction, m_state_offsets, displacement_block); const mfem::Vector gravity_gradient_direction = - make_read_only_value_view(direction, m_state_true_offsets, gravity_gradient_block); + make_read_only_value_view(direction, m_state_offsets, gravity_gradient_block); const mfem::Vector gravity_potential_direction = - make_read_only_value_view(direction, m_state_true_offsets, gravity_potential_block); + make_read_only_value_view(direction, m_state_offsets, gravity_potential_block); + + const field::FieldDofMap &displacement_map = m_linearization_context.GetDisplacementMap(); + const field::FieldDofMap &flux_map = m_linearization_context.GetGravityGradientMap(); + const field::FieldDofMap &potential_map = m_linearization_context.GetGravityPotentialMap(); + + mfem::Vector displacement_direction_true(displacement_map.full_size()); + mfem::Vector gravity_gradient_direction_true(flux_map.full_size()); + mfem::Vector gravity_potential_direction_true(potential_map.full_size()); + displacement_map.scatter(displacement_direction, displacement_direction_true); + flux_map.scatter(gravity_gradient_direction, gravity_gradient_direction_true); + potential_map.scatter(gravity_potential_direction, gravity_potential_direction_true); action.SetSize(Height()); action = 0.0; mfem::Vector gravity_gradient_action = - make_residual_view(action, m_residual_true_offsets, gravity_gradient_residual_block); + make_residual_view(action, m_residual_offsets, gravity_gradient_residual_block); mfem::Vector gravity_poisson_action = - make_residual_view(action, m_residual_true_offsets, gravity_poisson_residual_block); + make_residual_view(action, m_residual_offsets, gravity_poisson_residual_block); - mfem::Vector transpose_divergence_action; + mfem::Vector transpose_divergence_action_true; + mfem::Vector transpose_divergence_action(flux_map.reduced_size()); + mfem::Vector divergence_action_true; mfem::Vector source_action; - mfem::Vector mass_variation_action; - mfem::Vector source_variation_action; + mfem::Vector mass_variation_action_true; + mfem::Vector mass_variation_action(flux_map.reduced_size()); + mfem::Vector source_variation_action_true; + mfem::Vector source_variation_action(potential_map.reduced_size()); geometry_context.GetMassOperator().Mult(gravity_gradient_direction, gravity_gradient_action); geometry_context.GetSourceOperator().Mult(density_direction, source_action); kernels::apply_mapped_hdiv_mass_variation( - m_fem, m_domain_mapper, gravity_gradient, displacement, displacement_direction, mass_variation_action + m_fem, m_domain_mapper, gravity_gradient, displacement, displacement_direction_true, + mass_variation_action_true ); + flux_map.gather(mass_variation_action_true, mass_variation_action); kernels::apply_mapped_source_variation( - m_fem, m_domain_mapper, density, displacement, displacement_direction, source_variation_action + m_fem, m_domain_mapper, density, displacement, displacement_direction_true, source_variation_action_true ); + potential_map.gather(source_variation_action_true, source_variation_action); - transpose_divergence_action.SetSize(gravity_gradient_action.Size()); - m_fem.gravityContext.BT->Mult(gravity_potential_direction, transpose_divergence_action); + transpose_divergence_action_true.SetSize(flux_map.full_size()); + m_fem.gravityContext.BT->Mult(gravity_potential_direction_true, transpose_divergence_action_true); + flux_map.gather(transpose_divergence_action_true, transpose_divergence_action); gravity_gradient_action += transpose_divergence_action; gravity_gradient_action += mass_variation_action; - m_fem.gravityContext.b_form->Mult(gravity_gradient_direction, gravity_poisson_action); + divergence_action_true.SetSize(potential_map.full_size()); + m_fem.gravityContext.b_form->Mult(gravity_gradient_direction_true, divergence_action_true); + potential_map.gather(divergence_action_true, gravity_poisson_action); gravity_poisson_action -= source_action; gravity_poisson_action -= source_variation_action; diff --git a/libmeanfield/impl/operators/prepared_displacement_operator.cpp b/libmeanfield/impl/operators/prepared_displacement_operator.cpp index da7dcbb..8e5078c 100644 --- a/libmeanfield/impl/operators/prepared_displacement_operator.cpp +++ b/libmeanfield/impl/operators/prepared_displacement_operator.cpp @@ -153,10 +153,11 @@ namespace mean_field::operators { ); } - const mfem::Vector &density = m_gravityContext.GetDensity(); - const mfem::Vector &displacement = m_gravityContext.GetGeometryContext().GetDisplacement(); + const mfem::Vector density = m_gravityContext.GetDensityMap().gather(m_gravityContext.GetDensityTrue()); + const mfem::Vector displacement = + m_gravityContext.GetDisplacementMap().gather(m_gravityContext.GetGeometryContext().GetDisplacementTrue()); - m_isPrepared = false; + m_isPrepared = false; PreparedDisplacementResidualReport report; @@ -170,13 +171,14 @@ namespace mean_field::operators { {.density = density, .displacement = displacement}, make_rotational_dependencies(dependencies), rotation ); - if (report.DidAnyChildWork() || m_cachedResidual.Size() != m_fem.displacementFes->GetTrueVSize()) { + if (report.DidAnyChildWork() || + m_cachedResidual.Size() != m_gravityContext.GetDisplacementMap().reduced_size()) { AssembleResidual(); report.assembledResidual = true; } MFEM_VERIFY( - m_cachedResidual.Size() == m_fem.displacementFes->GetTrueVSize(), + m_cachedResidual.Size() == m_gravityContext.GetDisplacementMap().reduced_size(), "PreparedDisplacementResidualOperator produced a cached " "residual with the wrong size." ); @@ -445,10 +447,13 @@ namespace mean_field::operators { utils::blocks::get_residual_block(utils::blocks::barotropic_constant_field.mass_normalization_term); MFEM_VERIFY( - m_layout.size(densityValue) == f.densityFes->GetTrueVSize() && - m_layout.size(displacementValue) == f.displacementFes->GetTrueVSize() && - m_layout.size(gravityGradientValue) == f.gravityFluxFes->GetTrueVSize() && - m_layout.size(gravityPotentialValue) == f.gravityPotentialFes->GetTrueVSize() && + m_layout.size(densityValue) == m_preparedOperator.GetGravityContext().GetDensityMap().reduced_size() && + m_layout.size(displacementValue) == + m_preparedOperator.GetGravityContext().GetDisplacementMap().reduced_size() && + m_layout.size(gravityGradientValue) == + m_preparedOperator.GetGravityContext().GetGravityGradientMap().reduced_size() && + m_layout.size(gravityPotentialValue) == + m_preparedOperator.GetGravityContext().GetGravityPotentialMap().reduced_size() && m_layout.size(barotropicConstantValue) == 1, "Prepared displacement-residual MFEM adapter received " "incompatible barotropic value-block sizes." @@ -461,11 +466,16 @@ namespace mean_field::operators { ); MFEM_VERIFY( - m_layout.size(gravityGradientResidual) == f.gravityFluxFes->GetTrueVSize() && - m_layout.size(gravityPotentialResidual) == f.gravityPotentialFes->GetTrueVSize() && - m_layout.size(densityResidual) == f.densityFes->GetTrueVSize() && - m_layout.size(displacementResidual) == f.displacementFes->GetTrueVSize() && - m_layout.size(enthalpyResidual) == f.enthalpyFes->GetTrueVSize() && m_layout.size(massResidual) == 1, + m_layout.size(gravityGradientResidual) == + m_preparedOperator.GetGravityContext().GetGravityGradientMap().reduced_size() && + m_layout.size(gravityPotentialResidual) == + m_preparedOperator.GetGravityContext().GetGravityPotentialMap().reduced_size() && + m_layout.size(densityResidual) == + m_preparedOperator.GetGravityContext().GetDensityMap().reduced_size() && + m_layout.size(displacementResidual) == + m_preparedOperator.GetGravityContext().GetDisplacementMap().reduced_size() && + m_layout.size(enthalpyResidual) == m_preparedOperator.GetPressureOperator().GetEnthalpySize() && + m_layout.size(massResidual) == 1, "Prepared displacement-residual MFEM adapter received " "incompatible barotropic residual-block sizes." ); diff --git a/libmeanfield/impl/operators/prepared_gravity_displacement_force.cpp b/libmeanfield/impl/operators/prepared_gravity_displacement_force.cpp index 34a3f00..024921f 100644 --- a/libmeanfield/impl/operators/prepared_gravity_displacement_force.cpp +++ b/libmeanfield/impl/operators/prepared_gravity_displacement_force.cpp @@ -54,9 +54,11 @@ namespace mean_field::operators { } kernels::apply_gravity_displacement_force_residual( - m_fem, m_domainMapper, m_gravityContext.GetDensity(), m_gravityContext.GetGravityGradient(), - m_gravityContext.GetGeometryContext().GetDisplacement(), m_cachedResidual + m_fem, m_domainMapper, m_gravityContext.GetDensityTrue(), m_gravityContext.GetGravityGradientTrue(), + m_gravityContext.GetGeometryContext().GetDisplacementTrue(), m_actionTrue ); + m_cachedResidual.SetSize(m_gravityContext.GetDisplacementMap().reduced_size()); + m_gravityContext.GetDisplacementMap().gather(m_actionTrue, m_cachedResidual); m_preparedRevisions = requestedRevisions; ++m_residualPreparationCount; @@ -77,10 +79,15 @@ namespace mean_field::operators { ) const { VerifyPrepared(); + m_densityVariationTrue.SetSize(m_gravityContext.GetDensityMap().full_size()); + m_gravityContext.GetDensityMap().scatter(densityVariation, m_densityVariationTrue); + kernels::apply_gravity_displacement_force_density_action( - m_fem, m_domainMapper, densityVariation, m_gravityContext.GetGravityGradient(), - m_gravityContext.GetGeometryContext().GetDisplacement(), action + m_fem, m_domainMapper, m_densityVariationTrue, m_gravityContext.GetGravityGradientTrue(), + m_gravityContext.GetGeometryContext().GetDisplacementTrue(), m_actionTrue ); + action.SetSize(m_gravityContext.GetDisplacementMap().reduced_size()); + m_gravityContext.GetDisplacementMap().gather(m_actionTrue, action); ++m_densityJacobianStatistics.applications; } @@ -91,10 +98,15 @@ namespace mean_field::operators { ) const { VerifyPrepared(); + m_gravityGradientVariationTrue.SetSize(m_gravityContext.GetGravityGradientMap().full_size()); + m_gravityContext.GetGravityGradientMap().scatter(gravityGradientVariation, m_gravityGradientVariationTrue); + kernels::apply_gravity_displacement_force_gradient_action( - m_fem, m_domainMapper, m_gravityContext.GetDensity(), gravityGradientVariation, - m_gravityContext.GetGeometryContext().GetDisplacement(), action + m_fem, m_domainMapper, m_gravityContext.GetDensityTrue(), m_gravityGradientVariationTrue, + m_gravityContext.GetGeometryContext().GetDisplacementTrue(), m_actionTrue ); + action.SetSize(m_gravityContext.GetDisplacementMap().reduced_size()); + m_gravityContext.GetDisplacementMap().gather(m_actionTrue, action); ++m_gravityGradientJacobianStatistics.applications; } @@ -105,10 +117,15 @@ namespace mean_field::operators { ) const { VerifyPrepared(); + m_displacementVariationTrue.SetSize(m_gravityContext.GetDisplacementMap().full_size()); + m_gravityContext.GetDisplacementMap().scatter(displacementVariation, m_displacementVariationTrue); + kernels::apply_gravity_displacement_force_displacement_action( - m_fem, m_domainMapper, m_gravityContext.GetDensity(), m_gravityContext.GetGravityGradient(), - displacementVariation, m_gravityContext.GetGeometryContext().GetDisplacement(), action + m_fem, m_domainMapper, m_gravityContext.GetDensityTrue(), m_gravityContext.GetGravityGradientTrue(), + m_displacementVariationTrue, m_gravityContext.GetGeometryContext().GetDisplacementTrue(), m_actionTrue ); + action.SetSize(m_gravityContext.GetDisplacementMap().reduced_size()); + m_gravityContext.GetDisplacementMap().gather(m_actionTrue, action); ++m_displacementJacobianStatistics.applications; } @@ -121,11 +138,20 @@ namespace mean_field::operators { ) const { VerifyPrepared(); + m_densityVariationTrue.SetSize(m_gravityContext.GetDensityMap().full_size()); + m_gravityGradientVariationTrue.SetSize(m_gravityContext.GetGravityGradientMap().full_size()); + m_displacementVariationTrue.SetSize(m_gravityContext.GetDisplacementMap().full_size()); + m_gravityContext.GetDensityMap().scatter(densityVariation, m_densityVariationTrue); + m_gravityContext.GetGravityGradientMap().scatter(gravityGradientVariation, m_gravityGradientVariationTrue); + m_gravityContext.GetDisplacementMap().scatter(displacementVariation, m_displacementVariationTrue); + kernels::apply_gravity_displacement_force_complete_action( - m_fem, m_domainMapper, m_gravityContext.GetDensity(), densityVariation, - m_gravityContext.GetGravityGradient(), gravityGradientVariation, displacementVariation, - m_gravityContext.GetGeometryContext().GetDisplacement(), action + m_fem, m_domainMapper, m_gravityContext.GetDensityTrue(), m_densityVariationTrue, + m_gravityContext.GetGravityGradientTrue(), m_gravityGradientVariationTrue, m_displacementVariationTrue, + m_gravityContext.GetGeometryContext().GetDisplacementTrue(), m_actionTrue ); + action.SetSize(m_gravityContext.GetDisplacementMap().reduced_size()); + m_gravityContext.GetDisplacementMap().gather(m_actionTrue, action); ++m_densityJacobianStatistics.applications; ++m_gravityGradientJacobianStatistics.applications; @@ -196,8 +222,6 @@ namespace mean_field::operators { ), m_layout(layout), m_preparedOperator(preparedOperator) { - const fem::FEM &f = m_preparedOperator.GetFEM(); - using Form = utils::blocks::barotropic_equilibrium_form; constexpr auto densityValue = utils::blocks::get_value_block(utils::blocks::density_field.mass_term); @@ -212,10 +236,13 @@ namespace mean_field::operators { utils::blocks::get_residual_block(utils::blocks::displacement_field.geometry_term); MFEM_VERIFY( - m_layout.size(densityValue) == f.densityFes->GetTrueVSize() && - m_layout.size(displacementValue) == f.displacementFes->GetTrueVSize() && - m_layout.size(gravityGradientValue) == f.gravityFluxFes->GetTrueVSize() && - m_layout.size(displacementResidual) == f.displacementFes->GetTrueVSize(), + m_layout.size(densityValue) == m_preparedOperator.GetGravityContext().GetDensityMap().reduced_size() && + m_layout.size(displacementValue) == + m_preparedOperator.GetGravityContext().GetDisplacementMap().reduced_size() && + m_layout.size(gravityGradientValue) == + m_preparedOperator.GetGravityContext().GetGravityGradientMap().reduced_size() && + m_layout.size(displacementResidual) == + m_preparedOperator.GetGravityContext().GetDisplacementMap().reduced_size(), "Prepared gravity-displacement-force MFEM adapter received " "incompatible coupled block sizes." ); @@ -287,4 +314,4 @@ namespace mean_field::operators { const GravityDisplacementForceLayout &PreparedGravityDisplacementForceJacobianOperator::GetLayout() const noexcept { return m_layout; } -} // namespace mean_field::operators \ No newline at end of file +} // namespace mean_field::operators diff --git a/libmeanfield/impl/operators/prepared_gravity_source.cpp b/libmeanfield/impl/operators/prepared_gravity_source.cpp index e70251e..554b645 100644 --- a/libmeanfield/impl/operators/prepared_gravity_source.cpp +++ b/libmeanfield/impl/operators/prepared_gravity_source.cpp @@ -9,13 +9,16 @@ module mean_field; import :operators.prepared_gravity_source; namespace { + using DomainSchema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema; + int get_operator_height(const mean_field::fem::FEM &f) { MFEM_VERIFY( f.gravityPotentialFes != nullptr, "PreparedMappedGravitySourceOperator requires the " "gravity-potential " "finite-element space." ); - return f.gravityPotentialFes->GetTrueVSize(); + return mean_field::field::make_field_dof_map(*f.gravityPotentialFes) + .reduced_size(); } int get_operator_width(const mean_field::fem::FEM &f) { @@ -23,7 +26,8 @@ namespace { f.densityFes != nullptr, "PreparedMappedGravitySourceOperator requires the density " "finite-element space." ); - return f.densityFes->GetTrueVSize(); + return mean_field::field::make_field_dof_map(*f.densityFes) + .reduced_size(); } void true_to_local( @@ -238,7 +242,22 @@ namespace mean_field::operators { get_operator_width(f) ), m_fem(f), - m_domain_mapper(domain_mapper) { + m_domain_mapper(domain_mapper), + m_density_map( + field::make_field_dof_map< + field::Density, + DomainSchema>(*f.densityFes) + ), + m_potential_map( + field::make_field_dof_map< + field::Gravity, + DomainSchema>(*f.gravityPotentialFes) + ), + m_displacement_map( + field::make_field_dof_map< + field::Displacement, + DomainSchema>(*f.displacementFes) + ) { MFEM_VERIFY(f.mesh != nullptr, "PreparedMappedGravitySourceOperator requires a mesh."); MFEM_VERIFY( f.densityFes != nullptr, "PreparedMappedGravitySourceOperator requires the density " @@ -275,26 +294,28 @@ namespace mean_field::operators { utils::populate_element_mask(f.mesh.get(), utils::DOMAINS::STELLAR, m_stellar_marker); } - void PreparedMappedGravitySourceOperator::Prepare(const mfem::Vector &displacement_true) { + void PreparedMappedGravitySourceOperator::Prepare(const mfem::Vector &displacement) { MFEM_VERIFY( - displacement_true.Size() == m_fem.displacementFes->GetTrueVSize(), + displacement.Size() == m_displacement_map.reduced_size(), "PreparedMappedGravitySourceOperator received a displacement " "vector " "with the wrong size." ); - for (int i = 0; i < displacement_true.Size(); ++i) { + for (int i = 0; i < displacement.Size(); ++i) { MFEM_VERIFY( - std::isfinite(displacement_true(i)), "PreparedMappedGravitySourceOperator received a non-finite " - "displacement value." + std::isfinite(displacement(i)), "PreparedMappedGravitySourceOperator received a non-finite " + "displacement value." ); } m_is_prepared = false; + m_displacement_true.SetSize(m_displacement_map.full_size()); + m_displacement_map.scatter(displacement, m_displacement_true); m_elements.clear(); m_elements.reserve(m_fem.mesh->GetNE()); - FrozenMappedGravitySourceCoefficient source_coefficient(m_fem, m_domain_mapper, displacement_true); + FrozenMappedGravitySourceCoefficient source_coefficient(m_fem, m_domain_mapper, m_displacement_true); for (int element_id = 0; element_id < m_fem.mesh->GetNE(); ++element_id) { const int attribute = m_fem.mesh->GetAttribute(element_id); @@ -379,7 +400,7 @@ namespace mean_field::operators { ++m_preparation_count; } void PreparedMappedGravitySourceOperator::Mult( - const mfem::Vector &density_true, + const mfem::Vector &density, mfem::Vector &action ) const { MFEM_VERIFY( @@ -388,13 +409,16 @@ namespace mean_field::operators { ); MFEM_VERIFY( - density_true.Size() == Width(), "PreparedMappedGravitySourceOperator received a density vector " - "with the wrong size." + density.Size() == Width(), "PreparedMappedGravitySourceOperator received a density vector " + "with the wrong size." ); + m_density_true.SetSize(m_density_map.full_size()); + m_density_map.scatter(density, m_density_true); + mfem::Vector density_local; - true_to_local(*m_fem.densityFes, density_true, density_local); + true_to_local(*m_fem.densityFes, m_density_true, density_local); mfem::Vector local_action(m_fem.gravityPotentialFes->GetVSize()); local_action = 0.0; @@ -432,11 +456,13 @@ namespace mean_field::operators { local_action.AddElementVector(data.potential_dofs, element_action); } - local_to_true(*m_fem.gravityPotentialFes, local_action, action); + local_to_true(*m_fem.gravityPotentialFes, local_action, m_action_true); + action.SetSize(Height()); + m_potential_map.gather(m_action_true, action); } void PreparedMappedGravitySourceOperator::MultTranspose( - const mfem::Vector &potential_true, + const mfem::Vector &potential, mfem::Vector &action ) const { MFEM_VERIFY( @@ -445,13 +471,16 @@ namespace mean_field::operators { ); MFEM_VERIFY( - potential_true.Size() == Height(), "PreparedMappedGravitySourceOperator received a potential vector " - "with the wrong size." + potential.Size() == Height(), "PreparedMappedGravitySourceOperator received a potential vector " + "with the wrong size." ); + m_potential_true.SetSize(m_potential_map.full_size()); + m_potential_map.scatter(potential, m_potential_true); + mfem::Vector potential_local; - true_to_local(*m_fem.gravityPotentialFes, potential_true, potential_local); + true_to_local(*m_fem.gravityPotentialFes, m_potential_true, potential_local); mfem::Vector local_action(m_fem.densityFes->GetVSize()); local_action = 0.0; @@ -486,7 +515,9 @@ namespace mean_field::operators { local_action.AddElementVector(data.density_dofs, element_action); } - local_to_true(*m_fem.densityFes, local_action, action); + local_to_true(*m_fem.densityFes, local_action, m_action_true); + action.SetSize(Width()); + m_density_map.gather(m_action_true, action); } bool PreparedMappedGravitySourceOperator::IsPrepared() const noexcept { return m_is_prepared; @@ -495,4 +526,16 @@ namespace mean_field::operators { std::uint64_t PreparedMappedGravitySourceOperator::GetPreparationCount() const noexcept { return m_preparation_count; } + + const field::FieldDofMap &PreparedMappedGravitySourceOperator::GetDensityMap() const noexcept { + return m_density_map; + } + + const field::FieldDofMap &PreparedMappedGravitySourceOperator::GetPotentialMap() const noexcept { + return m_potential_map; + } + + const field::FieldDofMap &PreparedMappedGravitySourceOperator::GetDisplacementMap() const noexcept { + return m_displacement_map; + } } // namespace mean_field::operators diff --git a/libmeanfield/impl/operators/prepared_hdiv_mass.cpp b/libmeanfield/impl/operators/prepared_hdiv_mass.cpp index f469990..e84e022 100644 --- a/libmeanfield/impl/operators/prepared_hdiv_mass.cpp +++ b/libmeanfield/impl/operators/prepared_hdiv_mass.cpp @@ -8,12 +8,15 @@ module mean_field; import :operators.prepared_hdiv_mass; namespace { + using DomainSchema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema; + int get_operator_size(const mean_field::fem::FEM &f) { MFEM_VERIFY( f.gravityFluxFes != nullptr, "PreparedMappedHDivMassOperator requires the " "gravity-gradient finite-element space." ); - return f.gravityFluxFes->GetTrueVSize(); + return mean_field::field::make_field_dof_map(*f.gravityFluxFes) + .reduced_size(); } void true_to_local( @@ -222,7 +225,17 @@ namespace mean_field::operators { ) : Operator(get_operator_size(f)), m_fem(f), - m_domain_mapper(domain_mapper) { + m_domain_mapper(domain_mapper), + m_flux_map( + field::make_field_dof_map< + field::Gravity, + DomainSchema>(*f.gravityFluxFes) + ), + m_displacement_map( + field::make_field_dof_map< + field::Displacement, + DomainSchema>(*f.displacementFes) + ) { MFEM_VERIFY(f.mesh != nullptr, "PreparedMappedHDivMassOperator requires a mesh."); MFEM_VERIFY( f.gravityFluxFes != nullptr, "PreparedMappedHDivMassOperator requires the " @@ -269,22 +282,25 @@ namespace mean_field::operators { validate_uniform_domain_discretization(f, m_vacuum_marker, vacuum_element_id); } - void PreparedMappedHDivMassOperator::Prepare(const mfem::Vector &displacement_true) { + void PreparedMappedHDivMassOperator::Prepare(const mfem::Vector &displacement) { MFEM_VERIFY( - displacement_true.Size() == m_fem.displacementFes->GetTrueVSize(), + displacement.Size() == m_displacement_map.reduced_size(), "PreparedMappedHDivMassOperator received a displacement vector " "with " "the wrong size." ); - for (int i = 0; i < displacement_true.Size(); ++i) { + for (int i = 0; i < displacement.Size(); ++i) { MFEM_VERIFY( - std::isfinite(displacement_true(i)), "PreparedMappedHDivMassOperator received a non-finite " - "displacement " - "value." + std::isfinite(displacement(i)), "PreparedMappedHDivMassOperator received a non-finite " + "displacement " + "value." ); } + m_displacement_true.SetSize(m_displacement_map.full_size()); + m_displacement_map.scatter(displacement, m_displacement_true); + const int stellar_element_id = find_representative_element(m_fem, m_stellar_marker); const int vacuum_element_id = find_representative_element(m_fem, m_vacuum_marker); @@ -299,9 +315,9 @@ namespace mean_field::operators { m_vacuum_mass_coefficient.reset(); m_stellar_mass_coefficient = - std::make_unique(m_fem, m_domain_mapper, displacement_true, false); + std::make_unique(m_fem, m_domain_mapper, m_displacement_true, false); m_vacuum_mass_coefficient = - std::make_unique(m_fem, m_domain_mapper, displacement_true, true); + std::make_unique(m_fem, m_domain_mapper, m_displacement_true, true); m_mass_form = std::make_unique(m_fem.gravityFluxFes.get()); m_mass_form->SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL); @@ -328,7 +344,7 @@ namespace mean_field::operators { } void PreparedMappedHDivMassOperator::Mult( - const mfem::Vector &gravity_gradient_true, + const mfem::Vector &gravity_gradient, mfem::Vector &action ) const { MFEM_VERIFY( @@ -340,13 +356,16 @@ namespace mean_field::operators { "assembled partial-assembly form." ); MFEM_VERIFY( - gravity_gradient_true.Size() == Width(), - "PreparedMappedHDivMassOperator received a gravity-gradient vector " - "with the wrong size." + gravity_gradient.Size() == Width(), "PreparedMappedHDivMassOperator received a gravity-gradient vector " + "with the wrong size." ); + m_flux_true.SetSize(m_flux_map.full_size()); + m_action_true.SetSize(m_flux_map.full_size()); + m_flux_map.scatter(gravity_gradient, m_flux_true); + m_mass_form->Mult(m_flux_true, m_action_true); action.SetSize(Height()); - m_mass_form->Mult(gravity_gradient_true, action); + m_flux_map.gather(m_action_true, action); } bool PreparedMappedHDivMassOperator::IsPrepared() const noexcept { @@ -356,4 +375,12 @@ namespace mean_field::operators { std::uint64_t PreparedMappedHDivMassOperator::GetPreparationCount() const noexcept { return m_preparation_count; } -} // namespace mean_field::operators \ No newline at end of file + + const field::FieldDofMap &PreparedMappedHDivMassOperator::GetFluxMap() const noexcept { + return m_flux_map; + } + + const field::FieldDofMap &PreparedMappedHDivMassOperator::GetDisplacementMap() const noexcept { + return m_displacement_map; + } +} // namespace mean_field::operators diff --git a/libmeanfield/impl/operators/prepared_hydrostatic_equilibrium.cpp b/libmeanfield/impl/operators/prepared_hydrostatic_equilibrium.cpp index 7f4f6fb..cc76503 100644 --- a/libmeanfield/impl/operators/prepared_hydrostatic_equilibrium.cpp +++ b/libmeanfield/impl/operators/prepared_hydrostatic_equilibrium.cpp @@ -12,6 +12,8 @@ module mean_field; import :operators.prepared_hydrostatic_equilibrium; namespace { + using DomainSchema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema; + void true_to_local( const mfem::ParFiniteElementSpace &finiteElementSpace, const mfem::Vector &trueVector, @@ -151,9 +153,18 @@ namespace mean_field::operators { "displacement finite-element space." ); - m_enthalpySize = f.enthalpyFes->GetTrueVSize(); - m_gravityPotentialSize = f.gravityPotentialFes->GetTrueVSize(); - m_displacementSize = f.displacementFes->GetTrueVSize(); + const field::FieldDofMap enthalpyMap = + field::make_field_dof_map(*f.enthalpyFes); + + const field::FieldDofMap gravityPotentialMap = + field::make_field_dof_map(*f.gravityPotentialFes); + + const field::FieldDofMap displacementMap = + field::make_field_dof_map(*f.displacementFes); + + m_enthalpySize = enthalpyMap.reduced_size(); + m_gravityPotentialSize = gravityPotentialMap.reduced_size(); + m_displacementSize = displacementMap.reduced_size(); m_residualSize = m_enthalpySize; m_totalSize = m_enthalpySize + m_gravityPotentialSize + 1 + m_displacementSize; @@ -265,6 +276,11 @@ namespace mean_field::operators { m_domainMapper.GetDimension() == m_fem.mesh->Dimension(), "The hydrostatic operator's stateless mapper " "dimension does not match the mesh dimension." ); + + m_enthalpyVariationTrue.SetSize(m_context.GetEnthalpyMap().full_size()); + m_gravityPotentialVariationTrue.SetSize(m_context.GetGravityPotentialMap().full_size()); + m_displacementVariationTrue.SetSize(m_context.GetDisplacementMap().full_size()); + m_fullEnthalpyAction.SetSize(m_context.GetEnthalpyMap().full_size()); } PreparedHydrostaticEquilibriumReport PreparedHydrostaticEquilibriumOperator::Prepare( @@ -320,8 +336,8 @@ namespace mean_field::operators { ); MFEM_VERIFY( - m_cachedResidual.Size() == m_fem.enthalpyFes->GetTrueVSize(), - "The prepared hydrostatic residual has the wrong size." + m_cachedResidual.Size() == m_context.GetEnthalpyMap().reduced_size(), + "The prepared hydrostatic residual has the wrong supported size." ); m_isPrepared = true; @@ -708,7 +724,10 @@ namespace mean_field::operators { localResidual.AddElementVector(data.enthalpyDofs, elementResidual); } - local_to_true(*m_fem.enthalpyFes, localResidual, m_cachedResidual); + local_to_true(*m_fem.enthalpyFes, localResidual, m_fullEnthalpyAction); + + m_cachedResidual.SetSize(m_context.GetEnthalpyMap().reduced_size()); + m_context.GetEnthalpyMap().gather(m_fullEnthalpyAction, m_cachedResidual); } void PreparedHydrostaticEquilibriumOperator::BuildResidual(mfem::Vector &residual) const { @@ -723,9 +742,16 @@ namespace mean_field::operators { ) const { VerifyPrepared(); + MFEM_VERIFY( + enthalpyVariation.Size() == m_context.GetEnthalpyMap().reduced_size(), + "Prepared hydrostatic enthalpy variation has the wrong supported size." + ); + + m_context.GetEnthalpyMap().scatter(enthalpyVariation, m_enthalpyVariationTrue); + mfem::Vector enthalpyVariationLocal; - true_to_local(*m_fem.enthalpyFes, enthalpyVariation, enthalpyVariationLocal); + true_to_local(*m_fem.enthalpyFes, m_enthalpyVariationTrue, enthalpyVariationLocal); mfem::Vector localAction(m_fem.enthalpyFes->GetVSize()); @@ -752,7 +778,10 @@ namespace mean_field::operators { localAction.AddElementVector(data.enthalpyDofs, elementAction); } - local_to_true(*m_fem.enthalpyFes, localAction, action); + local_to_true(*m_fem.enthalpyFes, localAction, m_fullEnthalpyAction); + + action.SetSize(m_context.GetEnthalpyMap().reduced_size()); + m_context.GetEnthalpyMap().gather(m_fullEnthalpyAction, action); ++m_algebraicJacobianStatistics.enthalpyApplications; } @@ -763,9 +792,18 @@ namespace mean_field::operators { ) const { VerifyPrepared(); + MFEM_VERIFY( + gravityPotentialVariation.Size() == m_context.GetGravityPotentialMap().reduced_size(), + "Prepared hydrostatic gravity-potential variation has the wrong supported size." + ); + + m_context.GetGravityPotentialMap().scatter(gravityPotentialVariation, m_gravityPotentialVariationTrue); + mfem::Vector gravityPotentialVariationLocal; - true_to_local(*m_fem.gravityPotentialFes, gravityPotentialVariation, gravityPotentialVariationLocal); + true_to_local( + *m_fem.gravityPotentialFes, m_gravityPotentialVariationTrue, gravityPotentialVariationLocal + ); mfem::Vector localAction(m_fem.enthalpyFes->GetVSize()); @@ -792,7 +830,10 @@ namespace mean_field::operators { localAction.AddElementVector(data.enthalpyDofs, elementAction); } - local_to_true(*m_fem.enthalpyFes, localAction, action); + local_to_true(*m_fem.enthalpyFes, localAction, m_fullEnthalpyAction); + + action.SetSize(m_context.GetEnthalpyMap().reduced_size()); + m_context.GetEnthalpyMap().gather(m_fullEnthalpyAction, action); ++m_algebraicJacobianStatistics.gravityPotentialApplications; } @@ -824,7 +865,10 @@ namespace mean_field::operators { localAction.AddElementVector(data.enthalpyDofs, elementAction); } - local_to_true(*m_fem.enthalpyFes, localAction, action); + local_to_true(*m_fem.enthalpyFes, localAction, m_fullEnthalpyAction); + + action.SetSize(m_context.GetEnthalpyMap().reduced_size()); + m_context.GetEnthalpyMap().gather(m_fullEnthalpyAction, action); ++m_algebraicJacobianStatistics.bernoulliConstantApplications; } @@ -842,12 +886,27 @@ namespace mean_field::operators { "a non-finite Bernoulli-constant variation." ); + MFEM_VERIFY( + enthalpyVariation.Size() == m_context.GetEnthalpyMap().reduced_size(), + "Prepared hydrostatic algebraic enthalpy variation has the wrong supported size." + ); + + MFEM_VERIFY( + gravityPotentialVariation.Size() == m_context.GetGravityPotentialMap().reduced_size(), + "Prepared hydrostatic algebraic gravity-potential variation has the wrong supported size." + ); + + m_context.GetEnthalpyMap().scatter(enthalpyVariation, m_enthalpyVariationTrue); + m_context.GetGravityPotentialMap().scatter(gravityPotentialVariation, m_gravityPotentialVariationTrue); + mfem::Vector enthalpyVariationLocal; mfem::Vector gravityPotentialVariationLocal; - true_to_local(*m_fem.enthalpyFes, enthalpyVariation, enthalpyVariationLocal); + true_to_local(*m_fem.enthalpyFes, m_enthalpyVariationTrue, enthalpyVariationLocal); - true_to_local(*m_fem.gravityPotentialFes, gravityPotentialVariation, gravityPotentialVariationLocal); + true_to_local( + *m_fem.gravityPotentialFes, m_gravityPotentialVariationTrue, gravityPotentialVariationLocal + ); mfem::Vector localAction(m_fem.enthalpyFes->GetVSize()); @@ -898,7 +957,10 @@ namespace mean_field::operators { localAction.AddElementVector(data.enthalpyDofs, elementAction); } - local_to_true(*m_fem.enthalpyFes, localAction, action); + local_to_true(*m_fem.enthalpyFes, localAction, m_fullEnthalpyAction); + + action.SetSize(m_context.GetEnthalpyMap().reduced_size()); + m_context.GetEnthalpyMap().gather(m_fullEnthalpyAction, action); ++m_algebraicJacobianStatistics.combinedApplications; } @@ -909,9 +971,16 @@ namespace mean_field::operators { ) const { VerifyPrepared(); + MFEM_VERIFY( + displacementVariation.Size() == m_context.GetDisplacementMap().reduced_size(), + "Prepared hydrostatic displacement variation has the wrong supported size." + ); + + m_context.GetDisplacementMap().scatter(displacementVariation, m_displacementVariationTrue); + mfem::Vector displacementVariationLocal; - true_to_local(*m_fem.displacementFes, displacementVariation, displacementVariationLocal); + true_to_local(*m_fem.displacementFes, m_displacementVariationTrue, displacementVariationLocal); mfem::Vector localAction(m_fem.enthalpyFes->GetVSize()); @@ -1018,7 +1087,10 @@ namespace mean_field::operators { localAction.AddElementVector(data.enthalpyDofs, elementAction); } - local_to_true(*m_fem.enthalpyFes, localAction, action); + local_to_true(*m_fem.enthalpyFes, localAction, m_fullEnthalpyAction); + + action.SetSize(m_context.GetEnthalpyMap().reduced_size()); + m_context.GetEnthalpyMap().gather(m_fullEnthalpyAction, action); ++m_displacementJacobianStatistics.applications; } @@ -1087,6 +1159,18 @@ namespace mean_field::operators { return m_fem; } + const field::FieldDofMap &PreparedHydrostaticEquilibriumOperator::GetEnthalpyMap() const noexcept { + return m_context.GetEnthalpyMap(); + } + + const field::FieldDofMap &PreparedHydrostaticEquilibriumOperator::GetGravityPotentialMap() const noexcept { + return m_context.GetGravityPotentialMap(); + } + + const field::FieldDofMap &PreparedHydrostaticEquilibriumOperator::GetDisplacementMap() const noexcept { + return m_context.GetDisplacementMap(); + } + void PreparedHydrostaticEquilibriumOperator::VerifyPrepared() const { MFEM_VERIFY( m_isPrepared, "PreparedHydrostaticEquilibriumOperator must be " @@ -1099,7 +1183,7 @@ namespace mean_field::operators { const PreparedHydrostaticEquilibriumOperator &preparedOperator ) : mfem::Operator( - f.enthalpyFes != nullptr ? f.enthalpyFes->GetTrueVSize() : 0, + HydrostaticJacobianBlockLayout(f).GetResidualSize(), HydrostaticJacobianBlockLayout(f).GetTotalSize() ), m_layout(f), diff --git a/libmeanfield/impl/operators/prepared_mass_normalization.cpp b/libmeanfield/impl/operators/prepared_mass_normalization.cpp index 856ec7e..3097225 100644 --- a/libmeanfield/impl/operators/prepared_mass_normalization.cpp +++ b/libmeanfield/impl/operators/prepared_mass_normalization.cpp @@ -114,6 +114,15 @@ namespace mean_field::operators { "PreparedMassNormalizationOperator received a mapper with the " "wrong dimension." ); + + MFEM_VERIFY( + m_gravityContext.GetDensityMap().full_size() == m_fem.densityFes->GetTrueVSize() && + m_gravityContext.GetDisplacementMap().full_size() == m_fem.displacementFes->GetTrueVSize(), + "PreparedMassNormalizationOperator received incompatible shared FieldDof maps." + ); + + m_densityVariationTrue.SetSize(m_gravityContext.GetDensityMap().full_size()); + m_displacementVariationTrue.SetSize(m_gravityContext.GetDisplacementMap().full_size()); } PreparedMassNormalizationReport PreparedMassNormalizationOperator::Prepare( @@ -167,12 +176,12 @@ namespace mean_field::operators { } if (refreshGeometry) { - RefreshGeometry(m_gravityContext.GetGeometryContext().GetDisplacement()); + RefreshGeometry(m_gravityContext.GetGeometryContext().GetDisplacementTrue()); report.refreshedGeometry = true; } if (refreshDensity) { - RefreshDensity(m_gravityContext.GetDensity()); + RefreshDensity(m_gravityContext.GetDensityTrue()); report.refreshedDensity = true; } @@ -476,8 +485,16 @@ namespace mean_field::operators { mfem::Vector &action ) const { VerifyPrepared(); + + MFEM_VERIFY( + densityVariation.Size() == m_gravityContext.GetDensityMap().reduced_size(), + "Mass-normalization density action received a supported vector with the wrong size." + ); + validate_finite_vector(densityVariation, "Mass-normalization density action received a non-finite value."); + m_gravityContext.GetDensityMap().scatter(densityVariation, m_densityVariationTrue); + action.SetSize(1); - action(0) = GlobalSum(EvaluateDensityActionLocal(densityVariation)); + action(0) = GlobalSum(EvaluateDensityActionLocal(m_densityVariationTrue)); ++m_actionStatistics.densityApplications; } @@ -486,8 +503,18 @@ namespace mean_field::operators { mfem::Vector &action ) const { VerifyPrepared(); + + MFEM_VERIFY( + displacementVariation.Size() == m_gravityContext.GetDisplacementMap().reduced_size(), + "Mass-normalization displacement action received a supported vector with the wrong size." + ); + validate_finite_vector( + displacementVariation, "Mass-normalization displacement action received a non-finite value." + ); + m_gravityContext.GetDisplacementMap().scatter(displacementVariation, m_displacementVariationTrue); + action.SetSize(1); - action(0) = GlobalSum(EvaluateDisplacementActionLocal(displacementVariation)); + action(0) = GlobalSum(EvaluateDisplacementActionLocal(m_displacementVariationTrue)); ++m_actionStatistics.displacementApplications; } @@ -498,8 +525,25 @@ namespace mean_field::operators { ) const { VerifyPrepared(); + MFEM_VERIFY( + densityVariation.Size() == m_gravityContext.GetDensityMap().reduced_size(), + "Mass-normalization complete action received a supported density vector with the wrong size." + ); + MFEM_VERIFY( + displacementVariation.Size() == m_gravityContext.GetDisplacementMap().reduced_size(), + "Mass-normalization complete action received a supported displacement vector with the wrong size." + ); + validate_finite_vector(densityVariation, "Mass-normalization complete action received a non-finite density."); + validate_finite_vector( + displacementVariation, "Mass-normalization complete action received a non-finite displacement." + ); + + m_gravityContext.GetDensityMap().scatter(densityVariation, m_densityVariationTrue); + m_gravityContext.GetDisplacementMap().scatter(displacementVariation, m_displacementVariationTrue); + const double localAction = - EvaluateDensityActionLocal(densityVariation) + EvaluateDisplacementActionLocal(displacementVariation); + EvaluateDensityActionLocal(m_densityVariationTrue) + + EvaluateDisplacementActionLocal(m_displacementVariationTrue); action.SetSize(1); action(0) = GlobalSum(localAction); @@ -607,18 +651,25 @@ namespace mean_field::operators { constexpr auto massResidual = utils::blocks::get_residual_block(utils::blocks::barotropic_constant_field.mass_normalization_term); + using DomainSchema = utils::domain::CoreEnvelopeVacuumDomainSchema; + + const auto &gravityContext = m_preparedOperator.GetGravityContext(); + + const field::FieldDofMap enthalpyMap = + field::make_field_dof_map(*f.enthalpyFes); + MFEM_VERIFY( - m_layout.size(densityValue) == f.densityFes->GetTrueVSize() && - m_layout.size(displacementValue) == f.displacementFes->GetTrueVSize() && - m_layout.size(gravityGradientValue) == f.gravityFluxFes->GetTrueVSize() && - m_layout.size(gravityPotentialValue) == f.gravityPotentialFes->GetTrueVSize() && - m_layout.size(enthalpyValue) == f.enthalpyFes->GetTrueVSize() && + m_layout.size(densityValue) == gravityContext.GetDensityMap().reduced_size() && + m_layout.size(displacementValue) == gravityContext.GetDisplacementMap().reduced_size() && + m_layout.size(gravityGradientValue) == gravityContext.GetGravityGradientMap().reduced_size() && + m_layout.size(gravityPotentialValue) == gravityContext.GetGravityPotentialMap().reduced_size() && + m_layout.size(enthalpyValue) == enthalpyMap.reduced_size() && m_layout.size(barotropicConstantValue) == 1 && - m_layout.size(gravityGradientResidual) == f.gravityFluxFes->GetTrueVSize() && - m_layout.size(gravityPotentialResidual) == f.gravityPotentialFes->GetTrueVSize() && - m_layout.size(densityResidual) == f.densityFes->GetTrueVSize() && - m_layout.size(displacementResidual) == f.displacementFes->GetTrueVSize() && - m_layout.size(enthalpyResidual) == f.enthalpyFes->GetTrueVSize() && m_layout.size(massResidual) == 1, + m_layout.size(gravityGradientResidual) == gravityContext.GetGravityGradientMap().reduced_size() && + m_layout.size(gravityPotentialResidual) == gravityContext.GetGravityPotentialMap().reduced_size() && + m_layout.size(densityResidual) == gravityContext.GetDensityMap().reduced_size() && + m_layout.size(displacementResidual) == gravityContext.GetDisplacementMap().reduced_size() && + m_layout.size(enthalpyResidual) == enthalpyMap.reduced_size() && m_layout.size(massResidual) == 1, "Prepared mass-normalization MFEM adapter received incompatible " "barotropic block sizes." ); @@ -665,4 +716,4 @@ namespace mean_field::operators { const MassNormalizationLayout &PreparedMassNormalizationJacobianOperator::GetLayout() const noexcept { return m_layout; } -} // namespace mean_field::operators \ No newline at end of file +} // namespace mean_field::operators diff --git a/libmeanfield/impl/operators/prepared_rotation_displacement_force.cpp b/libmeanfield/impl/operators/prepared_rotation_displacement_force.cpp index 23e06aa..1def16c 100644 --- a/libmeanfield/impl/operators/prepared_rotation_displacement_force.cpp +++ b/libmeanfield/impl/operators/prepared_rotation_displacement_force.cpp @@ -79,15 +79,17 @@ namespace mean_field::operators { if (report.contextReport.preparedBaseState) { kernels::apply_rotational_displacement_force_residual( m_fem, m_domainMapper, *m_rotation, m_context.GetBaseDensityTrue(), m_context.GetDisplacementTrue(), - m_cachedResidual + m_actionTrue ); + m_cachedResidual.SetSize(m_context.GetDisplacementMap().reduced_size()); + m_context.GetDisplacementMap().gather(m_actionTrue, m_cachedResidual); ++m_residualPreparationCount; report.preparedResidual = true; } MFEM_VERIFY( - m_cachedResidual.Size() == m_fem.displacementFes->GetTrueVSize(), + m_cachedResidual.Size() == m_context.GetDisplacementMap().reduced_size(), "The prepared rotational-displacement-force residual has the " "wrong size." ); @@ -110,9 +112,14 @@ namespace mean_field::operators { ) const { VerifyPrepared(); + m_densityVariationTrue.SetSize(m_context.GetDensityMap().full_size()); + m_context.GetDensityMap().scatter(densityVariation, m_densityVariationTrue); + kernels::apply_rotational_displacement_force_density_action( - m_fem, m_domainMapper, *m_rotation, densityVariation, m_context.GetDisplacementTrue(), action + m_fem, m_domainMapper, *m_rotation, m_densityVariationTrue, m_context.GetDisplacementTrue(), m_actionTrue ); + action.SetSize(m_context.GetDisplacementMap().reduced_size()); + m_context.GetDisplacementMap().gather(m_actionTrue, action); ++m_densityJacobianStatistics.applications; } @@ -123,10 +130,15 @@ namespace mean_field::operators { ) const { VerifyPrepared(); + m_displacementVariationTrue.SetSize(m_context.GetDisplacementMap().full_size()); + m_context.GetDisplacementMap().scatter(displacementVariation, m_displacementVariationTrue); + kernels::apply_rotational_displacement_force_displacement_action( - m_fem, m_domainMapper, *m_rotation, m_context.GetBaseDensityTrue(), displacementVariation, - m_context.GetDisplacementTrue(), action + m_fem, m_domainMapper, *m_rotation, m_context.GetBaseDensityTrue(), m_displacementVariationTrue, + m_context.GetDisplacementTrue(), m_actionTrue ); + action.SetSize(m_context.GetDisplacementMap().reduced_size()); + m_context.GetDisplacementMap().gather(m_actionTrue, action); ++m_displacementJacobianStatistics.applications; } @@ -138,10 +150,17 @@ namespace mean_field::operators { ) const { VerifyPrepared(); + m_densityVariationTrue.SetSize(m_context.GetDensityMap().full_size()); + m_displacementVariationTrue.SetSize(m_context.GetDisplacementMap().full_size()); + m_context.GetDensityMap().scatter(densityVariation, m_densityVariationTrue); + m_context.GetDisplacementMap().scatter(displacementVariation, m_displacementVariationTrue); + kernels::apply_rotational_displacement_force_complete_action( - m_fem, m_domainMapper, *m_rotation, m_context.GetBaseDensityTrue(), densityVariation, displacementVariation, - m_context.GetDisplacementTrue(), action + m_fem, m_domainMapper, *m_rotation, m_context.GetBaseDensityTrue(), m_densityVariationTrue, + m_displacementVariationTrue, m_context.GetDisplacementTrue(), m_actionTrue ); + action.SetSize(m_context.GetDisplacementMap().reduced_size()); + m_context.GetDisplacementMap().gather(m_actionTrue, action); ++m_densityJacobianStatistics.applications; ++m_displacementJacobianStatistics.applications; @@ -207,8 +226,6 @@ namespace mean_field::operators { ), m_layout(layout), m_preparedOperator(preparedOperator) { - const fem::FEM &f = m_preparedOperator.GetFEM(); - using Form = utils::blocks::barotropic_equilibrium_form; constexpr auto densityValue = utils::blocks::get_value_block(utils::blocks::density_field.mass_term); @@ -220,9 +237,11 @@ namespace mean_field::operators { utils::blocks::get_residual_block(utils::blocks::displacement_field.geometry_term); MFEM_VERIFY( - m_layout.size(densityValue) == f.densityFes->GetTrueVSize() && - m_layout.size(displacementValue) == f.displacementFes->GetTrueVSize() && - m_layout.size(displacementResidual) == f.displacementFes->GetTrueVSize(), + m_layout.size(densityValue) == m_preparedOperator.GetContext().GetDensityMap().reduced_size() && + m_layout.size(displacementValue) == + m_preparedOperator.GetContext().GetDisplacementMap().reduced_size() && + m_layout.size(displacementResidual) == + m_preparedOperator.GetContext().GetDisplacementMap().reduced_size(), "Prepared rotational-displacement-force MFEM adapter received " "incompatible coupled block sizes." ); diff --git a/libmeanfield/impl/operators/prepared_stellar_equilibrium.cpp b/libmeanfield/impl/operators/prepared_stellar_equilibrium.cpp index 323a23c..318a488 100644 --- a/libmeanfield/impl/operators/prepared_stellar_equilibrium.cpp +++ b/libmeanfield/impl/operators/prepared_stellar_equilibrium.cpp @@ -55,21 +55,29 @@ namespace { return {valueSizes, residualSizes}; } - [[nodiscard]] mfem::Array make_gravity_state_offsets(const mean_field::fem::FEM &f) { + [[nodiscard]] mfem::Array make_gravity_state_offsets( + const mean_field::field::FieldDofMap &densityMap, + const mean_field::field::FieldDofMap &displacementMap, + const mean_field::field::FieldDofMap &gravityFluxMap, + const mean_field::field::FieldDofMap &gravityPotentialMap + ) { mfem::Array offsets(5); offsets[0] = 0; - offsets[1] = offsets[0] + f.densityFes->GetTrueVSize(); - offsets[2] = offsets[1] + f.displacementFes->GetTrueVSize(); - offsets[3] = offsets[2] + f.gravityFluxFes->GetTrueVSize(); - offsets[4] = offsets[3] + f.gravityPotentialFes->GetTrueVSize(); + offsets[1] = offsets[0] + densityMap.reduced_size(); + offsets[2] = offsets[1] + displacementMap.reduced_size(); + offsets[3] = offsets[2] + gravityFluxMap.reduced_size(); + offsets[4] = offsets[3] + gravityPotentialMap.reduced_size(); return offsets; } - [[nodiscard]] mfem::Array make_gravity_residual_offsets(const mean_field::fem::FEM &f) { + [[nodiscard]] mfem::Array make_gravity_residual_offsets( + const mean_field::field::FieldDofMap &gravityFluxMap, + const mean_field::field::FieldDofMap &gravityPotentialMap + ) { mfem::Array offsets(3); offsets[0] = 0; - offsets[1] = f.gravityFluxFes->GetTrueVSize(); - offsets[2] = offsets[1] + f.gravityPotentialFes->GetTrueVSize(); + offsets[1] = gravityFluxMap.reduced_size(); + offsets[2] = offsets[1] + gravityPotentialMap.reduced_size(); return offsets; } @@ -295,8 +303,16 @@ namespace mean_field::operators { gravityPotentialMap, enthalpyMap )), - gravityStateOffsets(make_gravity_state_offsets(f)), - gravityResidualOffsets(make_gravity_residual_offsets(f)) { + gravityStateOffsets(make_gravity_state_offsets( + densityMap, + displacementMap, + gravityFluxMap, + gravityPotentialMap + )), + gravityResidualOffsets(make_gravity_residual_offsets( + gravityFluxMap, + gravityPotentialMap + )) { } }; @@ -386,12 +402,7 @@ namespace mean_field::operators { domainMapper, m_gravityContext ), - m_targetMass(targetMass), - m_densityMap(std::move(constructionData.densityMap)), - m_displacementMap(std::move(constructionData.displacementMap)), - m_gravityFluxMap(std::move(constructionData.gravityFluxMap)), - m_gravityPotentialMap(std::move(constructionData.gravityPotentialMap)), - m_enthalpyMap(std::move(constructionData.enthalpyMap)) { + m_targetMass(targetMass) { MFEM_VERIFY( std::isfinite(m_targetMass) && m_targetMass > 0.0, "PreparedStellarEquilibriumOperator requires a finite, positive target mass." @@ -402,35 +413,12 @@ namespace mean_field::operators { "PreparedStellarEquilibriumOperator has inconsistent block dimensions." ); - MFEM_VERIFY( - m_displacementMap.is_identity(), "PreparedStellarEquilibriumOperator currently requires Displacement " - "support to span the full MFEM true-DOF space." - ); - MFEM_VERIFY( - m_gravityFluxMap.is_identity(), "PreparedStellarEquilibriumOperator currently requires gravity-flux " - "support to span the full MFEM true-DOF space." - ); - MFEM_VERIFY( - m_gravityPotentialMap.is_identity(), "PreparedStellarEquilibriumOperator currently requires " - "gravity-potential support to span the full MFEM true-DOF space." - ); + m_gravityState.SetSize(m_gravityStateOffsets.Last()); - m_fullDensity.SetSize(m_densityMap.full_size()); - m_fullEnthalpy.SetSize(m_enthalpyMap.full_size()); - m_fullGravityState.SetSize(m_gravityStateOffsets.Last()); + m_gravityDirection.SetSize(m_gravityStateOffsets.Last()); - m_fullDensityVariation.SetSize(m_densityMap.full_size()); - m_fullEnthalpyVariation.SetSize(m_enthalpyMap.full_size()); - m_fullGravityDirection.SetSize(m_gravityStateOffsets.Last()); - m_fullEnthalpyAction.SetSize(m_enthalpyMap.full_size()); - - m_fullDensity = 0.0; - m_fullEnthalpy = 0.0; - m_fullGravityState = 0.0; - m_fullDensityVariation = 0.0; - m_fullEnthalpyVariation = 0.0; - m_fullGravityDirection = 0.0; - m_fullEnthalpyAction = 0.0; + m_gravityState = 0.0; + m_gravityDirection = 0.0; } PreparedStellarEquilibriumReport PreparedStellarEquilibriumOperator::Prepare( @@ -502,16 +490,13 @@ namespace mean_field::operators { const mfem::Vector reducedEnthalpy = make_value_view(state, m_layout, enthalpyValue); const mfem::Vector bernoulli = make_value_view(state, m_layout, bernoulliValue); - m_densityMap.scatter(reducedDensity, m_fullDensity); - m_enthalpyMap.scatter(reducedEnthalpy, m_fullEnthalpy); - pack_gravity_vector( - m_fullGravityState, m_gravityStateOffsets, m_fullDensity, displacement, gravityGradient, gravityPotential + m_gravityState, m_gravityStateOffsets, reducedDensity, displacement, gravityGradient, gravityPotential ); PreparedStellarEquilibriumReport report; - report.gravity = m_gravityOperator.Prepare(m_fullGravityState, make_gravity_revisions(dependencies)); + report.gravity = m_gravityOperator.Prepare(m_gravityState, make_gravity_revisions(dependencies)); report.barotropicClosure = m_barotropicClosureOperator.Prepare( {.density = reducedDensity, .enthalpy = reducedEnthalpy, .displacement = displacement}, @@ -519,7 +504,7 @@ namespace mean_field::operators { ); report.hydrostatic = m_hydrostaticOperator.Prepare( - {.enthalpy = m_fullEnthalpy, + {.enthalpy = reducedEnthalpy, .gravityPotential = gravityPotential, .displacement = displacement, .bernoulliConstant = bernoulli(0)}, @@ -563,12 +548,13 @@ namespace mean_field::operators { mfem::Vector gravity; mfem::Vector closure; mfem::Vector displacement; + mfem::Vector hydrostatic; mfem::Vector mass; - m_gravityOperator.Mult(m_fullGravityState, gravity); + m_gravityOperator.Mult(m_gravityState, gravity); m_barotropicClosureOperator.BuildResidual(closure); m_displacementOperator.BuildResidual(displacement); - m_hydrostaticOperator.BuildResidual(m_fullEnthalpyAction); + m_hydrostaticOperator.BuildResidual(hydrostatic); m_massNormalizationOperator.BuildResidual(mass); m_cachedResidual.SetSize(Height()); @@ -600,10 +586,10 @@ namespace mean_field::operators { "The displacement residual has the wrong size." ); - { - mfem::Vector reducedEnthalpyResidual = make_residual_view(m_cachedResidual, m_layout, enthalpyResidual); - m_enthalpyMap.gather(m_fullEnthalpyAction, reducedEnthalpyResidual); - } + assign_residual_block( + m_cachedResidual, m_layout, enthalpyResidual, hydrostatic, + "The hydrostatic residual has the wrong size." + ); assign_residual_block( m_cachedResidual, m_layout, massResidual, mass, "The mass-normalization residual has the wrong size." @@ -665,37 +651,35 @@ namespace mean_field::operators { const mfem::Vector reducedEnthalpyDirection = make_value_view(direction, m_layout, enthalpyValue); const mfem::Vector bernoulliDirection = make_value_view(direction, m_layout, bernoulliValue); - m_densityMap.scatter(reducedDensityDirection, m_fullDensityVariation); - m_enthalpyMap.scatter(reducedEnthalpyDirection, m_fullEnthalpyVariation); - pack_gravity_vector( - m_fullGravityDirection, m_gravityStateOffsets, m_fullDensityVariation, displacementDirection, + m_gravityDirection, m_gravityStateOffsets, reducedDensityDirection, displacementDirection, gravityGradientDirection, gravityPotentialDirection ); mfem::Vector gravityAction; mfem::Vector closureAction; mfem::Vector displacementAction; + mfem::Vector hydrostaticAction; mfem::Vector massAction; - m_gravityJacobianOperator.Mult(m_fullGravityDirection, gravityAction); + m_gravityJacobianOperator.Mult(m_gravityDirection, gravityAction); m_barotropicClosureOperator.Mult( reducedDensityDirection, reducedEnthalpyDirection, displacementDirection, closureAction ); m_displacementOperator.ApplyCompleteJacobianAction( - m_fullDensityVariation, displacementDirection, gravityGradientDirection, reducedEnthalpyDirection, + reducedDensityDirection, displacementDirection, gravityGradientDirection, reducedEnthalpyDirection, displacementAction ); m_hydrostaticOperator.ApplyCompleteJacobianAction( - m_fullEnthalpyVariation, gravityPotentialDirection, bernoulliDirection(0), displacementDirection, - m_fullEnthalpyAction + reducedEnthalpyDirection, gravityPotentialDirection, bernoulliDirection(0), displacementDirection, + hydrostaticAction ); m_massNormalizationOperator.ApplyCompleteJacobianAction( - m_fullDensityVariation, displacementDirection, massAction + reducedDensityDirection, displacementDirection, massAction ); action.SetSize(Height()); @@ -727,10 +711,10 @@ namespace mean_field::operators { "The displacement Jacobian action has the wrong size." ); - { - mfem::Vector reducedEnthalpyAction = make_residual_view(action, m_layout, enthalpyResidual); - m_enthalpyMap.gather(m_fullEnthalpyAction, reducedEnthalpyAction); - } + assign_residual_block( + action, m_layout, enthalpyResidual, hydrostaticAction, + "The hydrostatic Jacobian action has the wrong size." + ); assign_residual_block( action, m_layout, massResidual, massAction, "The mass-normalization Jacobian action has the wrong size." diff --git a/libmeanfield/impl/physics/gravity.cpp b/libmeanfield/impl/physics/gravity.cpp index 1320975..ae5dc04 100644 --- a/libmeanfield/impl/physics/gravity.cpp +++ b/libmeanfield/impl/physics/gravity.cpp @@ -331,7 +331,10 @@ namespace mean_field::physics { auto source_form = std::make_unique(f, *f.domainMapperStateless); - source_form->Prepare(displacement_true); + using DomainSchema = utils::domain::CoreEnvelopeVacuumDomainSchema; + const field::FieldDofMap displacement_map = + field::make_field_dof_map(*f.displacementFes); + source_form->Prepare(displacement_map.gather(displacement_true)); f.gravityContext.source_form = std::move(source_form); // ========================================== @@ -440,13 +443,22 @@ namespace mean_field::physics { constexpr auto gravity_poisson_residual_block = utils::blocks::get_residual_block(utils::blocks::gravity_field.poisson_term); + using DomainSchema = utils::domain::CoreEnvelopeVacuumDomainSchema; + const field::FieldDofMap density_map = field::make_field_dof_map(*f.densityFes); + const field::FieldDofMap displacement_map = + field::make_field_dof_map(*f.displacementFes); + const field::FieldDofMap gravity_flux_map = + field::make_field_dof_map(*f.gravityFluxFes); + const field::FieldDofMap gravity_potential_map = + field::make_field_dof_map(*f.gravityPotentialFes); + const std::array value_sizes{ - f.densityFes->GetTrueVSize(), f.displacementFes->GetTrueVSize(), f.gravityFluxFes->GetTrueVSize(), - f.gravityPotentialFes->GetTrueVSize() + density_map.reduced_size(), displacement_map.reduced_size(), gravity_flux_map.reduced_size(), + gravity_potential_map.reduced_size() }; const std::array residual_sizes{ - f.gravityFluxFes->GetTrueVSize(), f.gravityPotentialFes->GetTrueVSize() + gravity_flux_map.reduced_size(), gravity_potential_map.reduced_size() }; const utils::blocks::form_layout layout(value_sizes, residual_sizes); @@ -457,6 +469,9 @@ namespace mean_field::physics { grid_function_to_true_dofs(*f.densityFes, rho, density_true); grid_function_to_true_dofs(*f.displacementFes, displacement, displacement_true); + const mfem::Vector density = density_map.gather(density_true); + const mfem::Vector reduced_displacement = displacement_map.gather(displacement_true); + operators::context::gravity_field::GravityFieldLinearizationContext linearization_context( f, *f.domainMapperStateless ); @@ -474,18 +489,18 @@ namespace mean_field::physics { ); operators::ReducedGravityFieldOperator reduced_operator( - gravity_operator, reduced_geometry_context, displacement_true + gravity_operator, reduced_geometry_context, reduced_displacement ); mfem::Vector right_hand_side; - reduced_operator.BuildRightHandSide(density_true, right_hand_side); + reduced_operator.BuildRightHandSide(density, right_hand_side); MFEM_VERIFY( right_hand_side.Size() == reduced_operator.Height(), "The reduced gravity right-hand side has the wrong size." ); - mfem::BlockVector gravity_state(reduced_operator.GetGravityTrueOffsets()); + mfem::BlockVector gravity_state(reduced_operator.GetGravityOffsets()); gravity_state = 0.0; mfem::MINRESSolver minres(f.mesh->GetComm()); @@ -501,9 +516,14 @@ namespace mean_field::physics { GravitySolution solution(f); - solution.gradPhi.SetFromTrueDofs(gravity_state.GetBlock(gravity_gradient_residual_block)); + const mfem::Vector gravity_flux_true = + gravity_flux_map.scatter(gravity_state.GetBlock(gravity_gradient_residual_block)); + const mfem::Vector gravity_potential_true = + gravity_potential_map.scatter(gravity_state.GetBlock(gravity_poisson_residual_block)); - solution.phi.SetFromTrueDofs(gravity_state.GetBlock(gravity_poisson_residual_block)); + solution.gradPhi.SetFromTrueDofs(gravity_flux_true); + + solution.phi.SetFromTrueDofs(gravity_potential_true); return solution; } diff --git a/libmeanfield/interface/operators/contexts/gravity_field_context.cppm b/libmeanfield/interface/operators/contexts/gravity_field_context.cppm index 62bef30..9ef3f03 100644 --- a/libmeanfield/interface/operators/contexts/gravity_field_context.cppm +++ b/libmeanfield/interface/operators/contexts/gravity_field_context.cppm @@ -6,6 +6,7 @@ module; export module mean_field:operators.context.gravity_field; export import :fem; +export import :field.mfem; export import :mapping.domain_mapper; export import :operators.prepared_gravity_source; export import :operators.prepared_hdiv_mass; @@ -69,14 +70,15 @@ export namespace mean_field::operators::context::gravity_field { GravityFieldGeometryContext &operator=(GravityFieldGeometryContext &&) = delete; GravityFieldGeometryPreparation Prepare( - const mfem::Vector &displacement_true, + const mfem::Vector &displacement, DiscretizationRevision discretization_revision, DisplacementRevision displacement_revision ); [[nodiscard]] const PreparedMappedHDivMassOperator &GetMassOperator() const; [[nodiscard]] const PreparedMappedGravitySourceOperator &GetSourceOperator() const; - [[nodiscard]] const mfem::Vector &GetDisplacement() const; + [[nodiscard]] const mfem::Vector &GetDisplacementTrue() const; + [[nodiscard]] const field::FieldDofMap &GetDisplacementMap() const noexcept; [[nodiscard]] DiscretizationRevision GetDiscretizationRevision() const noexcept; [[nodiscard]] DisplacementRevision GetDisplacementRevision() const noexcept; [[nodiscard]] bool IsPrepared() const noexcept; @@ -88,6 +90,7 @@ export namespace mean_field::operators::context::gravity_field { std::unique_ptr m_mass_operator; std::unique_ptr m_source_operator; + field::FieldDofMap m_displacement_map; mfem::Vector m_displacement_true; DiscretizationRevision m_discretization_revision; @@ -124,8 +127,12 @@ export namespace mean_field::operators::context::gravity_field { ); [[nodiscard]] const GravityFieldGeometryContext &GetGeometryContext() const; - [[nodiscard]] const mfem::Vector &GetDensity() const; - [[nodiscard]] const mfem::Vector &GetGravityGradient() const; + [[nodiscard]] const mfem::Vector &GetDensityTrue() const; + [[nodiscard]] const mfem::Vector &GetGravityGradientTrue() const; + [[nodiscard]] const field::FieldDofMap &GetDensityMap() const noexcept; + [[nodiscard]] const field::FieldDofMap &GetDisplacementMap() const noexcept; + [[nodiscard]] const field::FieldDofMap &GetGravityGradientMap() const noexcept; + [[nodiscard]] const field::FieldDofMap &GetGravityPotentialMap() const noexcept; [[nodiscard]] const GravityFieldRevisions &GetRevisions() const; [[nodiscard]] bool IsPrepared() const noexcept; @@ -134,10 +141,13 @@ export namespace mean_field::operators::context::gravity_field { GravityFieldGeometryContext m_geometry_context; + field::FieldDofMap m_density_map; + field::FieldDofMap m_gravity_gradient_map; + field::FieldDofMap m_gravity_potential_map; mfem::Vector m_density_true; mfem::Vector m_gravity_gradient_true; GravityFieldRevisions m_revisions; bool m_is_prepared{false}; }; -} // namespace mean_field::operators::context::gravity_field \ No newline at end of file +} // namespace mean_field::operators::context::gravity_field diff --git a/libmeanfield/interface/operators/contexts/hydrostatic_equilibrium_context.cppm b/libmeanfield/interface/operators/contexts/hydrostatic_equilibrium_context.cppm index 4c4596e..d5d8b26 100644 --- a/libmeanfield/interface/operators/contexts/hydrostatic_equilibrium_context.cppm +++ b/libmeanfield/interface/operators/contexts/hydrostatic_equilibrium_context.cppm @@ -8,6 +8,7 @@ module; export module mean_field:operators.context.hydrostatic_equilibrium; export import :fem; +export import :field.mfem; export import :mapping.domain_mapper; export namespace mean_field::operators::context::hydrostatic { @@ -52,6 +53,12 @@ export namespace mean_field::operators::context::hydrostatic { constexpr auto operator<=>(const HydrostaticEquilibriumDependencies &) const = default; }; + /* + * Frozen solver-facing state in registered FieldDof coordinates. + * + * The context owns the only expansion into MFEM true-DOF vectors used + * by the stateless hydrostatic kernels and prepared quadrature data. + */ struct HydrostaticEquilibriumStateView { const mfem::Vector &enthalpy; const mfem::Vector &gravityPotential; @@ -113,6 +120,12 @@ export namespace mean_field::operators::context::hydrostatic { [[nodiscard]] const HydrostaticPreparationStatistics &GetPreparationStatistics() const noexcept; + [[nodiscard]] const field::FieldDofMap &GetEnthalpyMap() const noexcept; + + [[nodiscard]] const field::FieldDofMap &GetGravityPotentialMap() const noexcept; + + [[nodiscard]] const field::FieldDofMap &GetDisplacementMap() const noexcept; + [[nodiscard]] const mfem::Vector &GetBaseEnthalpyTrue() const; [[nodiscard]] const mfem::Vector &GetBaseGravityPotentialTrue() const; @@ -127,6 +140,10 @@ export namespace mean_field::operators::context::hydrostatic { const fem::FEM &m_f; const mapping::DomainMapperStateless &m_domainMapper; + field::FieldDofMap m_enthalpyMap; + field::FieldDofMap m_gravityPotentialMap; + field::FieldDofMap m_displacementMap; + mfem::Vector m_baseEnthalpyTrue; mfem::Vector m_baseGravityPotentialTrue; mfem::Vector m_displacementTrue; diff --git a/libmeanfield/interface/operators/contexts/rotation_displacement_force_context.cppm b/libmeanfield/interface/operators/contexts/rotation_displacement_force_context.cppm index 85be92c..95971d6 100644 --- a/libmeanfield/interface/operators/contexts/rotation_displacement_force_context.cppm +++ b/libmeanfield/interface/operators/contexts/rotation_displacement_force_context.cppm @@ -8,6 +8,7 @@ module; export module mean_field:operators.context.rotational_displacement_force; export import :fem; +export import :field.mfem; export import :mapping.domain_mapper; export namespace mean_field::operators::context::rotational_displacement_force { @@ -108,12 +109,16 @@ export namespace mean_field::operators::context::rotational_displacement_force { [[nodiscard]] const mfem::Vector &GetBaseDensityTrue() const; [[nodiscard]] const mfem::Vector &GetDisplacementTrue() const; + [[nodiscard]] const field::FieldDofMap &GetDensityMap() const noexcept; + [[nodiscard]] const field::FieldDofMap &GetDisplacementMap() const noexcept; private: void VerifyPrepared() const; const fem::FEM &m_f; + field::FieldDofMap m_densityMap; + field::FieldDofMap m_displacementMap; mfem::Vector m_baseDensityTrue; mfem::Vector m_displacementTrue; diff --git a/libmeanfield/interface/operators/gravity_field.cppm b/libmeanfield/interface/operators/gravity_field.cppm index 98a668c..24f5173 100644 --- a/libmeanfield/interface/operators/gravity_field.cppm +++ b/libmeanfield/interface/operators/gravity_field.cppm @@ -1,5 +1,4 @@ module; -#include #include #include @@ -10,21 +9,13 @@ export import :operators.gravity_field_jacobian; export import :operators.context.gravity_field; export namespace mean_field::operators { - enum class GravityResidualBlock : std::uint8_t { gradient_equation = 0, poisson_equation = 1, count = 2 }; - - constexpr int gravity_residual_block_index(const GravityResidualBlock block) noexcept { - return static_cast(block); - } - - inline constexpr int gravity_residual_block_count = gravity_residual_block_index(GravityResidualBlock::count); - class GravityFieldOperator final : public mfem::Operator { public: GravityFieldOperator( fem::FEM &f, const mapping::DomainMapperStateless &domain_mapper, context::gravity_field::GravityFieldLinearizationContext &linearization_context, - const mfem::Array &state_true_offsets, + const mfem::Array &state_offsets, GravityFieldJacobianOperator &jacobian ); @@ -40,9 +31,9 @@ export namespace mean_field::operators { Operator &GetGradient(const mfem::Vector &state) const override; - [[nodiscard]] const mfem::Array &GetStateTrueOffsets() const noexcept; + [[nodiscard]] const mfem::Array &GetStateOffsets() const noexcept; - [[nodiscard]] const mfem::Array &GetResidualTrueOffsets() const noexcept; + [[nodiscard]] const mfem::Array &GetResidualOffsets() const noexcept; [[nodiscard]] context::gravity_field::GravityFieldLinearizationContext &GetLinearizationContext() noexcept; @@ -66,8 +57,8 @@ export namespace mean_field::operators { fem::FEM &m_fem; const mapping::DomainMapperStateless &m_domain_mapper; context::gravity_field::GravityFieldLinearizationContext &m_linearization_context; - mfem::Array m_state_true_offsets; - mfem::Array m_residual_true_offsets; + mfem::Array m_state_offsets; + mfem::Array m_residual_offsets; GravityFieldJacobianOperator &m_jacobian; }; @@ -106,7 +97,7 @@ export namespace mean_field::operators { [[nodiscard]] const context::gravity_field::GravityFieldGeometryContext &GetGeometryContext() const noexcept; - [[nodiscard]] const mfem::Array &GetGravityTrueOffsets() const noexcept; + [[nodiscard]] const mfem::Array &GetGravityOffsets() const noexcept; private: void ValidateDisplacement(const mfem::Vector &displacement) const; @@ -117,7 +108,8 @@ export namespace mean_field::operators { private: GravityFieldOperator &m_gravity_field_operator; - mfem::Array m_gravity_true_offsets; + mfem::Array m_gravity_offsets; context::gravity_field::GravityFieldGeometryContext &m_gravity_field_geometry_context; + mfem::Vector m_displacement; }; -} // namespace mean_field::operators \ No newline at end of file +} // namespace mean_field::operators diff --git a/libmeanfield/interface/operators/gravity_field_jacobian.cppm b/libmeanfield/interface/operators/gravity_field_jacobian.cppm index 44ed8a1..2a0f7b0 100644 --- a/libmeanfield/interface/operators/gravity_field_jacobian.cppm +++ b/libmeanfield/interface/operators/gravity_field_jacobian.cppm @@ -13,8 +13,8 @@ export namespace mean_field::operators { fem::FEM &f, const mapping::DomainMapperStateless &domain_mapper, const context::gravity_field::GravityFieldLinearizationContext &linearization_context, - const mfem::Array &state_true_offsets, - const mfem::Array &residual_true_offsets + const mfem::Array &state_offsets, + const mfem::Array &residual_offsets ); void Mult( @@ -29,7 +29,7 @@ export namespace mean_field::operators { fem::FEM &m_fem; const mapping::DomainMapperStateless &m_domain_mapper; const context::gravity_field::GravityFieldLinearizationContext &m_linearization_context; - mfem::Array m_state_true_offsets; - mfem::Array m_residual_true_offsets; + mfem::Array m_state_offsets; + mfem::Array m_residual_offsets; }; -} // namespace mean_field::operators \ No newline at end of file +} // namespace mean_field::operators diff --git a/libmeanfield/interface/operators/kernels/gravity_kernels.cppm b/libmeanfield/interface/operators/kernels/gravity_kernels.cppm index 538dcba..d5a1025 100644 --- a/libmeanfield/interface/operators/kernels/gravity_kernels.cppm +++ b/libmeanfield/interface/operators/kernels/gravity_kernels.cppm @@ -5,6 +5,13 @@ export import :mapping.domain_mapper; export import :fem; export namespace mean_field::operators::kernels { + /* + * Kernel boundary contract + * ------------------------ + * These low-level element kernels operate on MFEM true-DOF vectors. The + * prepared operators and gravity contexts own FieldDof maps and are the + * only layer that gathers/scatters reduced public coordinates. + */ void apply_mapped_hdiv_mass( const fem::FEM &f, diff --git a/libmeanfield/interface/operators/prepared_gravity_displacement_force.cppm b/libmeanfield/interface/operators/prepared_gravity_displacement_force.cppm index 8bc7e6f..90be52c 100644 --- a/libmeanfield/interface/operators/prepared_gravity_displacement_force.cppm +++ b/libmeanfield/interface/operators/prepared_gravity_displacement_force.cppm @@ -113,6 +113,11 @@ export namespace mean_field::operators { context::gravity_field::GravityFieldRevisions m_preparedRevisions; mfem::Vector m_cachedResidual; + mutable mfem::Vector m_densityVariationTrue; + mutable mfem::Vector m_gravityGradientVariationTrue; + mutable mfem::Vector m_displacementVariationTrue; + mutable mfem::Vector m_actionTrue; + std::uint64_t m_residualPreparationCount{0}; mutable std::uint64_t m_residualApplicationCount{0}; diff --git a/libmeanfield/interface/operators/prepared_gravity_source.cppm b/libmeanfield/interface/operators/prepared_gravity_source.cppm index 6cc29f0..9d6b9a8 100644 --- a/libmeanfield/interface/operators/prepared_gravity_source.cppm +++ b/libmeanfield/interface/operators/prepared_gravity_source.cppm @@ -6,6 +6,7 @@ module; export module mean_field:operators.prepared_gravity_source; export import :fem; +export import :field.mfem; export import :mapping.domain_mapper; export namespace mean_field::operators { @@ -16,17 +17,21 @@ export namespace mean_field::operators { const mapping::DomainMapperStateless &domain_mapper ); - void Prepare(const mfem::Vector &displacement_true); + void Prepare(const mfem::Vector &displacement); void Mult( - const mfem::Vector &density_true, + const mfem::Vector &density, mfem::Vector &action ) const override; [[nodiscard]] bool IsPrepared() const noexcept; [[nodiscard]] std::uint64_t GetPreparationCount() const noexcept; + [[nodiscard]] const field::FieldDofMap &GetDensityMap() const noexcept; + [[nodiscard]] const field::FieldDofMap &GetPotentialMap() const noexcept; + [[nodiscard]] const field::FieldDofMap &GetDisplacementMap() const noexcept; + void MultTranspose( - const mfem::Vector &potential_true, + const mfem::Vector &potential, mfem::Vector &action ) const override; @@ -52,9 +57,18 @@ export namespace mean_field::operators { const fem::FEM &m_fem; const mapping::DomainMapperStateless &m_domain_mapper; + field::FieldDofMap m_density_map; + field::FieldDofMap m_potential_map; + field::FieldDofMap m_displacement_map; + mfem::Array m_stellar_marker; std::vector m_elements; + mutable mfem::Vector m_density_true; + mutable mfem::Vector m_potential_true; + mutable mfem::Vector m_action_true; + mfem::Vector m_displacement_true; + std::uint64_t m_preparation_count{0}; bool m_is_prepared{false}; }; diff --git a/libmeanfield/interface/operators/prepared_hdiv_mass.cppm b/libmeanfield/interface/operators/prepared_hdiv_mass.cppm index e4f2285..00a489d 100644 --- a/libmeanfield/interface/operators/prepared_hdiv_mass.cppm +++ b/libmeanfield/interface/operators/prepared_hdiv_mass.cppm @@ -5,6 +5,7 @@ module; export module mean_field:operators.prepared_hdiv_mass; export import :fem; +export import :field.mfem; export import :mapping.domain_mapper; export namespace mean_field::operators { @@ -15,26 +16,35 @@ export namespace mean_field::operators { const mapping::DomainMapperStateless &domain_mapper ); - void Prepare(const mfem::Vector &displacement_true); + void Prepare(const mfem::Vector &displacement); void Mult( - const mfem::Vector &gravity_gradient_true, + const mfem::Vector &gravity_gradient, mfem::Vector &action ) const override; [[nodiscard]] bool IsPrepared() const noexcept; [[nodiscard]] std::uint64_t GetPreparationCount() const noexcept; + [[nodiscard]] const field::FieldDofMap &GetFluxMap() const noexcept; + [[nodiscard]] const field::FieldDofMap &GetDisplacementMap() const noexcept; + private: const fem::FEM &m_fem; const mapping::DomainMapperStateless &m_domain_mapper; + field::FieldDofMap m_flux_map; + field::FieldDofMap m_displacement_map; + mfem::Array m_stellar_marker; mfem::Array m_vacuum_marker; std::unique_ptr m_stellar_mass_coefficient; std::unique_ptr m_vacuum_mass_coefficient; std::unique_ptr m_mass_form; + mutable mfem::Vector m_flux_true; + mutable mfem::Vector m_action_true; + mfem::Vector m_displacement_true; std::uint64_t m_preparation_count{0}; bool m_is_prepared{false}; }; -} // namespace mean_field::operators \ No newline at end of file +} // namespace mean_field::operators diff --git a/libmeanfield/interface/operators/prepared_hydrostatic_equilibrium_operator.cppm b/libmeanfield/interface/operators/prepared_hydrostatic_equilibrium_operator.cppm index 2d86fc7..45ece03 100644 --- a/libmeanfield/interface/operators/prepared_hydrostatic_equilibrium_operator.cppm +++ b/libmeanfield/interface/operators/prepared_hydrostatic_equilibrium_operator.cppm @@ -160,6 +160,12 @@ export namespace mean_field::operators { [[nodiscard]] const fem::FEM &GetFEM() const noexcept; + [[nodiscard]] const field::FieldDofMap &GetEnthalpyMap() const noexcept; + + [[nodiscard]] const field::FieldDofMap &GetGravityPotentialMap() const noexcept; + + [[nodiscard]] const field::FieldDofMap &GetDisplacementMap() const noexcept; + private: struct ElementPAData { int elementId{-1}; @@ -219,6 +225,11 @@ export namespace mean_field::operators { std::vector m_elements; mfem::Vector m_cachedResidual; + mutable mfem::Vector m_enthalpyVariationTrue; + mutable mfem::Vector m_gravityPotentialVariationTrue; + mutable mfem::Vector m_displacementVariationTrue; + mutable mfem::Vector m_fullEnthalpyAction; + std::uint64_t m_residualPreparationCount{0}; mutable std::uint64_t m_residualApplicationCount{0}; diff --git a/libmeanfield/interface/operators/prepared_mass_normalization.cppm b/libmeanfield/interface/operators/prepared_mass_normalization.cppm index e32e539..4a8b4da 100644 --- a/libmeanfield/interface/operators/prepared_mass_normalization.cppm +++ b/libmeanfield/interface/operators/prepared_mass_normalization.cppm @@ -61,6 +61,8 @@ export namespace mean_field::operators { * Density and displacement are borrowed from the shared gravity-field * linearization context. This keeps the mass row on exactly the same * frozen state and geometry as the gravity and mechanical rows. + * Jacobian directions use the corresponding registered FieldDof + * coordinates; expansion to MFEM true-DOF vectors is private. */ class PreparedMassNormalizationOperator final { public: @@ -153,6 +155,9 @@ export namespace mean_field::operators { MassNormalizationDependencies m_preparedDependencies; mfem::Vector m_cachedResidual; + mutable mfem::Vector m_densityVariationTrue; + mutable mfem::Vector m_displacementVariationTrue; + double m_currentMass{0.0}; double m_targetMass{0.0}; diff --git a/libmeanfield/interface/operators/prepared_rotation_displacement_force.cppm b/libmeanfield/interface/operators/prepared_rotation_displacement_force.cppm index ae37603..76cd209 100644 --- a/libmeanfield/interface/operators/prepared_rotation_displacement_force.cppm +++ b/libmeanfield/interface/operators/prepared_rotation_displacement_force.cppm @@ -111,6 +111,9 @@ export namespace mean_field::operators { std::optional m_rotation; mfem::Vector m_cachedResidual; + mutable mfem::Vector m_densityVariationTrue; + mutable mfem::Vector m_displacementVariationTrue; + mutable mfem::Vector m_actionTrue; context::rotational_displacement_force::RotationalDisplacementForceDependencies m_preparedDependencies; diff --git a/libmeanfield/interface/operators/prepared_stellar_equilibrium.cppm b/libmeanfield/interface/operators/prepared_stellar_equilibrium.cppm index 0c67311..ab5dd9d 100644 --- a/libmeanfield/interface/operators/prepared_stellar_equilibrium.cppm +++ b/libmeanfield/interface/operators/prepared_stellar_equilibrium.cppm @@ -158,20 +158,8 @@ export namespace mean_field::operators { mutable PreparedStellarEquilibriumStatistics m_statistics; bool m_isPrepared{false}; - field::FieldDofMap m_densityMap; - field::FieldDofMap m_displacementMap; - field::FieldDofMap m_gravityFluxMap; - field::FieldDofMap m_gravityPotentialMap; - field::FieldDofMap m_enthalpyMap; + mfem::Vector m_gravityState; - mfem::Vector m_fullDensity; - mfem::Vector m_fullEnthalpy; - mfem::Vector m_fullGravityState; - - mutable mfem::Vector m_fullDensityVariation; - mutable mfem::Vector m_fullEnthalpyVariation; - mutable mfem::Vector m_fullGravityDirection; - - mutable mfem::Vector m_fullEnthalpyAction; + mutable mfem::Vector m_gravityDirection; }; } // namespace mean_field::operators diff --git a/tests/operators/contexts/gravity_field_context.cpp b/tests/operators/contexts/gravity_field_context.cpp index a7fa840..65f7d87 100644 --- a/tests/operators/contexts/gravity_field_context.cpp +++ b/tests/operators/contexts/gravity_field_context.cpp @@ -10,18 +10,19 @@ namespace gravity_context = operators::context::gravity_field; TEST_CASE( "Gravity Field Linearization Context Applies Selective Invalidation", - tags::integration &tags::gravity &tags::contexts + tags::gravity_context ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); gravity_context::GravityFieldLinearizationContext context(f, *f.domainMapperStateless); - mfem::Vector density = prepared_test::make_deterministic_vector(f.densityFes->GetTrueVSize(), 0.11); - mfem::Vector displacement = prepared_test::make_displacement(f, 0.0); - mfem::Vector gravity_gradient = prepared_test::make_deterministic_vector(f.gravityFluxFes->GetTrueVSize(), 0.37); + mfem::Vector density = prepared_test::make_deterministic_vector(context.GetDensityMap().reduced_size(), 0.11); + mfem::Vector displacement = context.GetDisplacementMap().gather(prepared_test::make_displacement(f, 0.0)); + mfem::Vector gravity_gradient = + prepared_test::make_deterministic_vector(context.GetGravityGradientMap().reduced_size(), 0.37); mfem::Vector gravity_potential = - prepared_test::make_deterministic_vector(f.gravityPotentialFes->GetTrueVSize(), 0.63); + prepared_test::make_deterministic_vector(context.GetGravityPotentialMap().reduced_size(), 0.63); gravity_context::GravityFieldRevisions revisions; @@ -72,7 +73,7 @@ TEST_CASE( CHECK(density_report.updated_density); CHECK_FALSE(density_report.updated_gravity_gradient); CHECK_FALSE(density_report.geometry.DidAnyWork()); - CHECK(context.GetDensity()(0) == density(0)); + CHECK(context.GetDensityTrue()(context.GetDensityMap().true_dof(0)) == density(0)); gravity_gradient(0) -= 0.4; ++revisions.gravity_gradient.value; @@ -82,9 +83,9 @@ TEST_CASE( CHECK_FALSE(gradient_report.updated_density); CHECK(gradient_report.updated_gravity_gradient); CHECK_FALSE(gradient_report.geometry.DidAnyWork()); - CHECK(context.GetGravityGradient()(0) == gravity_gradient(0)); + CHECK(context.GetGravityGradientTrue()(context.GetGravityGradientMap().true_dof(0)) == gravity_gradient(0)); - displacement = prepared_test::make_displacement(f, 1.0); + displacement = context.GetDisplacementMap().gather(prepared_test::make_displacement(f, 1.0)); ++revisions.displacement.value; const gravity_context::GravityFieldPreparationReport displacement_report = context.Prepare(make_state(), revisions); @@ -114,18 +115,19 @@ TEST_CASE( TEST_CASE( "Gravity Field Linearization Context Owns Frozen Base Fields", - tags::integration &tags::gravity &tags::contexts + tags::gravity_context ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); gravity_context::GravityFieldLinearizationContext context(f, *f.domainMapperStateless); - mfem::Vector density = prepared_test::make_deterministic_vector(f.densityFes->GetTrueVSize(), 0.13); - mfem::Vector displacement = prepared_test::make_displacement(f, 0.4); - mfem::Vector gravity_gradient = prepared_test::make_deterministic_vector(f.gravityFluxFes->GetTrueVSize(), 0.47); + mfem::Vector density = prepared_test::make_deterministic_vector(context.GetDensityMap().reduced_size(), 0.13); + mfem::Vector displacement = context.GetDisplacementMap().gather(prepared_test::make_displacement(f, 0.4)); + mfem::Vector gravity_gradient = + prepared_test::make_deterministic_vector(context.GetGravityGradientMap().reduced_size(), 0.47); mfem::Vector gravity_potential = - prepared_test::make_deterministic_vector(f.gravityPotentialFes->GetTrueVSize(), 0.71); + prepared_test::make_deterministic_vector(context.GetGravityPotentialMap().reduced_size(), 0.71); gravity_context::GravityFieldRevisions revisions; @@ -137,24 +139,24 @@ TEST_CASE( revisions ); - const mfem::Vector frozen_density = context.GetDensity(); - const mfem::Vector frozen_displacement = context.GetGeometryContext().GetDisplacement(); - const mfem::Vector frozen_gravity_gradient = context.GetGravityGradient(); + const mfem::Vector frozen_density = context.GetDensityTrue(); + const mfem::Vector frozen_displacement = context.GetGeometryContext().GetDisplacementTrue(); + const mfem::Vector frozen_gravity_gradient = context.GetGravityGradientTrue(); density = 0.0; displacement = 0.0; gravity_gradient = 0.0; gravity_potential = 0.0; - CHECK(prepared_test::relative_error(context.GetDensity(), frozen_density, f.mesh->GetComm()) == 0.0); + CHECK(prepared_test::relative_error(context.GetDensityTrue(), frozen_density, f.mesh->GetComm()) == 0.0); CHECK( prepared_test::relative_error( - context.GetGeometryContext().GetDisplacement(), frozen_displacement, f.displacementFes->GetComm() + context.GetGeometryContext().GetDisplacementTrue(), frozen_displacement, f.displacementFes->GetComm() ) == 0.0 ); CHECK( prepared_test::relative_error( - context.GetGravityGradient(), frozen_gravity_gradient, f.gravityFluxFes->GetComm() + context.GetGravityGradientTrue(), frozen_gravity_gradient, f.gravityFluxFes->GetComm() ) == 0.0 ); @@ -167,22 +169,22 @@ TEST_CASE( ); CHECK_FALSE(unchanged_revision_report.DidAnyWork()); - CHECK(prepared_test::relative_error(context.GetDensity(), frozen_density, f.mesh->GetComm()) == 0.0); + CHECK(prepared_test::relative_error(context.GetDensityTrue(), frozen_density, f.mesh->GetComm()) == 0.0); CHECK( prepared_test::relative_error( - context.GetGeometryContext().GetDisplacement(), frozen_displacement, f.displacementFes->GetComm() + context.GetGeometryContext().GetDisplacementTrue(), frozen_displacement, f.displacementFes->GetComm() ) == 0.0 ); CHECK( prepared_test::relative_error( - context.GetGravityGradient(), frozen_gravity_gradient, f.gravityFluxFes->GetComm() + context.GetGravityGradientTrue(), frozen_gravity_gradient, f.gravityFluxFes->GetComm() ) == 0.0 ); } TEST_CASE( "Gravity Field Geometry Contexts Have Independent Prepared State", - tags::integration &tags::gravity &tags::contexts + tags::gravity_context ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -190,10 +192,13 @@ TEST_CASE( gravity_context::GravityFieldGeometryContext first_context(f, *f.domainMapperStateless); gravity_context::GravityFieldGeometryContext second_context(f, *f.domainMapperStateless); - const mfem::Vector first_displacement = prepared_test::make_displacement(f, 0.0); - const mfem::Vector second_displacement = prepared_test::make_displacement(f, 1.0); - const mfem::Vector gravity_gradient = - prepared_test::make_deterministic_vector(f.gravityFluxFes->GetTrueVSize(), 0.35); + const mfem::Vector first_displacement = + first_context.GetDisplacementMap().gather(prepared_test::make_displacement(f, 0.0)); + const mfem::Vector second_displacement = + second_context.GetDisplacementMap().gather(prepared_test::make_displacement(f, 1.0)); + const mfem::Vector gravity_gradient = prepared_test::make_field_map(f).gather( + prepared_test::make_deterministic_vector(f.gravityFluxFes->GetTrueVSize(), 0.35) + ); first_context.Prepare(first_displacement, {.value = 0}, {.value = 0}); second_context.Prepare(second_displacement, {.value = 0}, {.value = 0}); @@ -205,7 +210,8 @@ TEST_CASE( first_context.GetMassOperator().Mult(gravity_gradient, first_action); second_context.GetMassOperator().Mult(gravity_gradient, second_action_before); - const mfem::Vector updated_first_displacement = prepared_test::make_displacement(f, 0.6); + const mfem::Vector updated_first_displacement = + first_context.GetDisplacementMap().gather(prepared_test::make_displacement(f, 0.6)); first_context.Prepare(updated_first_displacement, {.value = 0}, {.value = 1}); second_context.GetMassOperator().Mult(gravity_gradient, second_action_after); diff --git a/tests/operators/contexts/hydrostatic_equilibrium_context.cpp b/tests/operators/contexts/hydrostatic_equilibrium_context.cpp index 2c2cdec..8297909 100644 --- a/tests/operators/contexts/hydrostatic_equilibrium_context.cpp +++ b/tests/operators/contexts/hydrostatic_equilibrium_context.cpp @@ -40,7 +40,7 @@ namespace hydrostatic_context_test_utils { TEST_CASE( "Hydrostatic Context Applies Selective Invalidation", - tags::barotrope &tags::contexts &tags::hydro &tags::prepared &tags::unit + tags::barotrope_hydrostatic_context &tags::unit ) { auto args = test_utils::setup_args(); @@ -48,12 +48,15 @@ TEST_CASE( mean_field::operators::context::hydrostatic::HydrostaticEquilibriumContext context(f, *f.domainMapperStateless); - mfem::Vector enthalpy = gravity_prepared_test_utils::make_deterministic_vector(f.enthalpyFes->GetTrueVSize(), 0.17); + mfem::Vector enthalpy = + field_dof_test_utils::make_deterministic_supported_vector(*f.enthalpyFes, 0.17); mfem::Vector gravityPotential = - gravity_prepared_test_utils::make_deterministic_vector(f.gravityPotentialFes->GetTrueVSize(), 0.31); + field_dof_test_utils::make_deterministic_supported_vector( + *f.gravityPotentialFes, 0.31 + ); - mfem::Vector displacement = gravity_prepared_test_utils::make_displacement(f, 0.35); + mfem::Vector displacement = field_dof_test_utils::make_supported_displacement(f, 0.35); double bernoulliConstant = 0.73; @@ -145,7 +148,7 @@ TEST_CASE( CHECK_FALSE(enthalpyReport.updatedGravityPotential); CHECK_FALSE(enthalpyReport.updatedDisplacement); CHECK_FALSE(enthalpyReport.updatedBernoulliConstant); - CHECK(context.GetBaseEnthalpyTrue()(0) == enthalpy(0)); + CHECK(context.GetBaseEnthalpyTrue()(context.GetEnthalpyMap().true_dof(0)) == enthalpy(0)); ++dependencies.gravityPotential.revision; @@ -161,7 +164,9 @@ TEST_CASE( CHECK_FALSE(gravityPotentialReport.updatedDisplacement); CHECK_FALSE(gravityPotentialReport.updatedBernoulliConstant); - CHECK(context.GetBaseGravityPotentialTrue()(0) == gravityPotential(0)); + CHECK( + context.GetBaseGravityPotentialTrue()(context.GetGravityPotentialMap().true_dof(0)) == gravityPotential(0) + ); ++dependencies.bernoulliConstant.revision; @@ -212,7 +217,7 @@ TEST_CASE( CHECK(displacementReport.updatedDisplacement); CHECK_FALSE(displacementReport.updatedBernoulliConstant); - CHECK(context.GetDisplacementTrue()(0) == displacement(0)); + CHECK(context.GetDisplacementTrue()(context.GetDisplacementMap().true_dof(0)) == displacement(0)); ++dependencies.discretization.revision; @@ -240,7 +245,7 @@ TEST_CASE( TEST_CASE( "Hydrostatic Context Uses Identity In Every Dependency", - tags::barotrope &tags::contexts &tags::hydro &tags::prepared &tags::unit + tags::barotrope_hydrostatic_context &tags::unit ) { auto args = test_utils::setup_args(); @@ -248,12 +253,15 @@ TEST_CASE( mean_field::operators::context::hydrostatic::HydrostaticEquilibriumContext context(f, *f.domainMapperStateless); - mfem::Vector enthalpy = gravity_prepared_test_utils::make_deterministic_vector(f.enthalpyFes->GetTrueVSize(), 0.23); + mfem::Vector enthalpy = + field_dof_test_utils::make_deterministic_supported_vector(*f.enthalpyFes, 0.23); const mfem::Vector gravityPotential = - gravity_prepared_test_utils::make_deterministic_vector(f.gravityPotentialFes->GetTrueVSize(), 0.41); + field_dof_test_utils::make_deterministic_supported_vector( + *f.gravityPotentialFes, 0.41 + ); - const mfem::Vector displacement = gravity_prepared_test_utils::make_displacement(f, 0.60); + const mfem::Vector displacement = field_dof_test_utils::make_supported_displacement(f, 0.60); constexpr double bernoulliConstant = 0.81; @@ -288,7 +296,10 @@ TEST_CASE( ); CHECK_FALSE(sameStampReport.DidAnyWork()); - CHECK(context.GetBaseEnthalpyTrue()(0) == firstFrozenEnthalpy(0)); + CHECK( + context.GetBaseEnthalpyTrue()(context.GetEnthalpyMap().true_dof(0)) == + firstFrozenEnthalpy(context.GetEnthalpyMap().true_dof(0)) + ); ++dependencies.enthalpy.identity; dependencies.enthalpy.revision = 0; @@ -301,7 +312,7 @@ TEST_CASE( hydrostatic_context_test_utils::check_base_only(newIdentityReport); CHECK(newIdentityReport.updatedEnthalpy); - CHECK(context.GetBaseEnthalpyTrue()(0) == enthalpy(0)); + CHECK(context.GetBaseEnthalpyTrue()(context.GetEnthalpyMap().true_dof(0)) == enthalpy(0)); CHECK(context.MatchesDependencies(dependencies)); ++dependencies.rotation.identity; diff --git a/tests/operators/contexts/rotation_displacement_force_context.cpp b/tests/operators/contexts/rotation_displacement_force_context.cpp index b20c74d..fd0d949 100644 --- a/tests/operators/contexts/rotation_displacement_force_context.cpp +++ b/tests/operators/contexts/rotation_displacement_force_context.cpp @@ -53,7 +53,7 @@ namespace rotational_displacement_force_context_test_utils { TEST_CASE( "Rotational Displacement Force Context Applies Selective Invalidation", - tags::centrifugal &tags::contexts &tags::unit + tags::rotation_context_unit ) { mean_field::utils::Args args = test_utils::setup_args(); @@ -61,15 +61,16 @@ TEST_CASE( REQUIRE(f.okay()); - mfem::Vector density = rotational_displacement_force_context_test_utils::make_density(f, 0.83); - - mfem::Vector displacement = gravity_prepared_test_utils::make_displacement(f, 0.47); - - auto dependencies = rotational_displacement_force_context_test_utils::make_dependencies(); + auto dependencies = rotational_displacement_force_context_test_utils::make_dependencies(); rotational_displacement_force_context_test_utils::Context context(f, *f.domainMapperStateless); - const auto initialReport = context.Prepare({.density = density, .displacement = displacement}, dependencies); + mfem::Vector densityTrue = rotational_displacement_force_context_test_utils::make_density(f, 0.83); + mfem::Vector density = context.GetDensityMap().gather(densityTrue); + mfem::Vector displacementTrue = gravity_prepared_test_utils::make_displacement(f, 0.47); + mfem::Vector displacement = context.GetDisplacementMap().gather(displacementTrue); + + const auto initialReport = context.Prepare({.density = density, .displacement = displacement}, dependencies); REQUIRE(context.IsPrepared()); CHECK(context.MatchesDependencies(dependencies)); @@ -80,14 +81,18 @@ TEST_CASE( CHECK(initialReport.updatedDensity); CHECK(initialReport.updatedDisplacement); + const mfem::Vector expectedDensityTrue = context.GetDensityMap().scatter(density); + const mfem::Vector expectedDisplacementTrue = context.GetDisplacementMap().scatter(displacement); + CHECK( - rotational_displacement_force_context_test_utils::relative_difference(context.GetBaseDensityTrue(), density) < - 1.0e-14 + rotational_displacement_force_context_test_utils::relative_difference( + context.GetBaseDensityTrue(), expectedDensityTrue + ) < 1.0e-14 ); CHECK( rotational_displacement_force_context_test_utils::relative_difference( - context.GetDisplacementTrue(), displacement + context.GetDisplacementTrue(), expectedDisplacementTrue ) < 1.0e-14 ); @@ -95,7 +100,8 @@ TEST_CASE( CHECK_FALSE(unchangedReport.DidAnyWork()); - density = rotational_displacement_force_context_test_utils::make_density(f, 1.17); + densityTrue = rotational_displacement_force_context_test_utils::make_density(f, 1.17); + density = context.GetDensityMap().gather(densityTrue); ++dependencies.density.revision; @@ -110,7 +116,8 @@ TEST_CASE( const mfem::Vector densityAfterDensityRevision(context.GetBaseDensityTrue()); - displacement = gravity_prepared_test_utils::make_displacement(f, 0.81); + displacementTrue = gravity_prepared_test_utils::make_displacement(f, 0.81); + displacement = context.GetDisplacementMap().gather(displacementTrue); ++dependencies.displacement.revision; @@ -151,7 +158,7 @@ TEST_CASE( TEST_CASE( "Rotational Displacement Force Context Treats New Identities As New " "Dependency Streams", - tags::centrifugal &tags::contexts &tags::unit + tags::rotation_context_unit ) { mean_field::utils::Args args = test_utils::setup_args(); @@ -159,14 +166,15 @@ TEST_CASE( REQUIRE(f.okay()); - const mfem::Vector density = rotational_displacement_force_context_test_utils::make_density(f, 0.91); - - const mfem::Vector displacement = gravity_prepared_test_utils::make_displacement(f, 0.39); - - auto dependencies = rotational_displacement_force_context_test_utils::make_dependencies(); + auto dependencies = rotational_displacement_force_context_test_utils::make_dependencies(); rotational_displacement_force_context_test_utils::Context context(f, *f.domainMapperStateless); + const mfem::Vector density = + context.GetDensityMap().gather(rotational_displacement_force_context_test_utils::make_density(f, 0.91)); + const mfem::Vector displacement = + context.GetDisplacementMap().gather(gravity_prepared_test_utils::make_displacement(f, 0.39)); + context.Prepare({.density = density, .displacement = displacement}, dependencies); dependencies.density.identity += 1000; diff --git a/tests/operators/gravity_displacement_force.cpp b/tests/operators/gravity_displacement_force.cpp index 9f40ffb..5600c90 100644 --- a/tests/operators/gravity_displacement_force.cpp +++ b/tests/operators/gravity_displacement_force.cpp @@ -58,14 +58,25 @@ namespace gravity_displacement_force_test_utils { ); [[nodiscard]] mean_field::operators::GravityDisplacementForceLayout make_layout(const mean_field::fem::FEM &f) { + using DomainSchema = gravity_prepared_test_utils::DomainSchema; + + const auto densityMap = gravity_prepared_test_utils::make_field_map(f); + const auto displacementMap = gravity_prepared_test_utils::make_field_map(f); + const auto gravityFluxMap = + mean_field::field::make_field_dof_map(*f.gravityFluxFes); + const auto gravityPotentialMap = + mean_field::field::make_field_dof_map(*f.gravityPotentialFes); + const auto enthalpyMap = + mean_field::field::make_field_dof_map(*f.enthalpyFes); + const std::array valueSizes{ - f.densityFes->GetTrueVSize(), f.displacementFes->GetTrueVSize(), f.gravityFluxFes->GetTrueVSize(), - f.gravityPotentialFes->GetTrueVSize(), f.enthalpyFes->GetTrueVSize(), 1 + densityMap.reduced_size(), displacementMap.reduced_size(), gravityFluxMap.reduced_size(), + gravityPotentialMap.reduced_size(), enthalpyMap.reduced_size(), 1 }; const std::array residualSizes{ - f.gravityFluxFes->GetTrueVSize(), f.gravityPotentialFes->GetTrueVSize(), f.densityFes->GetTrueVSize(), - f.displacementFes->GetTrueVSize(), f.enthalpyFes->GetTrueVSize(), 1 + gravityFluxMap.reduced_size(), gravityPotentialMap.reduced_size(), densityMap.reduced_size(), + displacementMap.reduced_size(), enthalpyMap.reduced_size(), 1 }; return {valueSizes, residualSizes}; @@ -221,10 +232,10 @@ namespace gravity_displacement_force_test_utils { const mean_field::operators::context::gravity_field::GravityFieldRevisions &revisions ) { context.Prepare( - {.density = density, - .displacement = displacement, - .gravity_gradient = gravityGradient, - .gravity_potential = gravityPotential}, + {.density = context.GetDensityMap().gather(density), + .displacement = context.GetDisplacementMap().gather(displacement), + .gravity_gradient = context.GetGravityGradientMap().gather(gravityGradient), + .gravity_potential = context.GetGravityPotentialMap().gather(gravityPotential)}, revisions ); } @@ -314,7 +325,7 @@ namespace gravity_displacement_force_test_utils { TEST_CASE( "Gravity Displacement Force Query Includes Every Registered Operand", - tags::gravity &tags::quadrature &tags::unit + tags::gravity_unit ) { using DisplacementField = mean_field::field::Field; @@ -350,7 +361,7 @@ TEST_CASE( TEST_CASE( "Gravity Displacement Force Uses Positive Grad-Phi Sign And Excludes " "Vacuum", - tags::gravity &tags::integration &tags::accuracy + tags::gravity_kernel_accuracy ) { mean_field::utils::Args args = test_utils::setup_args(); @@ -413,7 +424,7 @@ TEST_CASE( TEST_CASE( "Prepared Gravity Displacement Force Reuses Shared Gravity Revisions", - tags::gravity &tags::prepared &tags::integration + tags::gravity_prepared ) { mean_field::utils::Args args = test_utils::setup_args(); @@ -457,9 +468,11 @@ TEST_CASE( f, *f.domainMapperStateless, density, gravityGradient, displacement, kernelResidual ); + const mfem::Vector kernelResidualReduced = gravityContext.GetDisplacementMap().gather(kernelResidual); + CHECK( gravity_displacement_force_test_utils::relative_difference( - preparedResidual, kernelResidual, f.mesh->GetComm() + preparedResidual, kernelResidualReduced, f.mesh->GetComm() ) < 2.0e-12 ); @@ -493,7 +506,7 @@ TEST_CASE( TEST_CASE( "Gravity Displacement Force Jacobian Matches All Columns And Centered " "Differences", - tags::gravity &tags::prepared &tags::jacobian &tags::accuracy + tags::gravity_prepared_jacobian_accuracy ) { mean_field::utils::Args args = test_utils::setup_args(); @@ -532,19 +545,24 @@ TEST_CASE( preparedOperator.Prepare(); + const mfem::Vector densityDirectionReduced = gravityContext.GetDensityMap().gather(densityDirection); + const mfem::Vector gravityGradientDirectionReduced = + gravityContext.GetGravityGradientMap().gather(gravityGradientDirection); + const mfem::Vector displacementDirectionReduced = gravityContext.GetDisplacementMap().gather(displacementDirection); + mfem::Vector densityAction; mfem::Vector gravityAction; mfem::Vector displacementAction; mfem::Vector completeAction; - preparedOperator.ApplyDensityJacobianAction(densityDirection, densityAction); + preparedOperator.ApplyDensityJacobianAction(densityDirectionReduced, densityAction); - preparedOperator.ApplyGravityGradientJacobianAction(gravityGradientDirection, gravityAction); + preparedOperator.ApplyGravityGradientJacobianAction(gravityGradientDirectionReduced, gravityAction); - preparedOperator.ApplyDisplacementJacobianAction(displacementDirection, displacementAction); + preparedOperator.ApplyDisplacementJacobianAction(displacementDirectionReduced, displacementAction); preparedOperator.ApplyCompleteJacobianAction( - densityDirection, displacementDirection, gravityGradientDirection, completeAction + densityDirectionReduced, displacementDirectionReduced, gravityGradientDirectionReduced, completeAction ); mfem::Vector summedColumns(densityAction); @@ -559,29 +577,34 @@ TEST_CASE( mfem::Vector zeroDensity(densityDirection.Size()); mfem::Vector zeroGravity(gravityGradientDirection.Size()); mfem::Vector zeroDisplacement(displacementDirection.Size()); - zeroDensity = 0.0; - zeroGravity = 0.0; - zeroDisplacement = 0.0; + zeroDensity = 0.0; + zeroGravity = 0.0; + zeroDisplacement = 0.0; - constexpr double step = 1.0e-5; + constexpr double step = 1.0e-5; - const mfem::Vector densityDifference = gravity_displacement_force_test_utils::centered_difference( + const mfem::Vector densityDifferenceTrue = gravity_displacement_force_test_utils::centered_difference( f, density, densityDirection, gravityGradient, zeroGravity, displacement, zeroDisplacement, step ); - const mfem::Vector gravityDifference = gravity_displacement_force_test_utils::centered_difference( + const mfem::Vector gravityDifferenceTrue = gravity_displacement_force_test_utils::centered_difference( f, density, zeroDensity, gravityGradient, gravityGradientDirection, displacement, zeroDisplacement, step ); - const mfem::Vector displacementDifference = gravity_displacement_force_test_utils::centered_difference( + const mfem::Vector displacementDifferenceTrue = gravity_displacement_force_test_utils::centered_difference( f, density, zeroDensity, gravityGradient, zeroGravity, displacement, displacementDirection, step ); - const mfem::Vector completeDifference = gravity_displacement_force_test_utils::centered_difference( + const mfem::Vector completeDifferenceTrue = gravity_displacement_force_test_utils::centered_difference( f, density, densityDirection, gravityGradient, gravityGradientDirection, displacement, displacementDirection, step ); + const mfem::Vector densityDifference = gravityContext.GetDisplacementMap().gather(densityDifferenceTrue); + const mfem::Vector gravityDifference = gravityContext.GetDisplacementMap().gather(gravityDifferenceTrue); + const mfem::Vector displacementDifference = gravityContext.GetDisplacementMap().gather(displacementDifferenceTrue); + const mfem::Vector completeDifference = gravityContext.GetDisplacementMap().gather(completeDifferenceTrue); + const double densityError = gravity_displacement_force_test_utils::relative_difference(densityAction, densityDifference, f.mesh->GetComm()); @@ -609,7 +632,7 @@ TEST_CASE( TEST_CASE( "Prepared Gravity Displacement Force MFEM Adapter Routes Only R-d", - tags::gravity &tags::prepared &tags::mfem_operators &tags::unit + tags::gravity_prepared_unit ) { mean_field::utils::Args args = test_utils::setup_args(); @@ -648,22 +671,27 @@ TEST_CASE( preparedOperator.Prepare(); - const auto layout = gravity_displacement_force_test_utils::make_layout(f); + const mfem::Vector densityDirectionReduced = gravityContext.GetDensityMap().gather(densityDirection); + const mfem::Vector gravityGradientDirectionReduced = + gravityContext.GetGravityGradientMap().gather(gravityGradientDirection); + const mfem::Vector displacementDirectionReduced = gravityContext.GetDisplacementMap().gather(displacementDirection); + + const auto layout = gravity_displacement_force_test_utils::make_layout(f); mean_field::operators::PreparedGravityDisplacementForceJacobianOperator adapter(layout, preparedOperator); mfem::BlockVector direction(layout.value_offsets()); - direction = 0.0; + direction = 0.0; - direction.GetBlock(gravity_displacement_force_test_utils::densityValue) = densityDirection; + direction.GetBlock(gravity_displacement_force_test_utils::densityValue) = densityDirectionReduced; - direction.GetBlock(gravity_displacement_force_test_utils::displacementValue) = displacementDirection; + direction.GetBlock(gravity_displacement_force_test_utils::displacementValue) = displacementDirectionReduced; - direction.GetBlock(gravity_displacement_force_test_utils::gravityGradientValue) = gravityGradientDirection; + direction.GetBlock(gravity_displacement_force_test_utils::gravityGradientValue) = gravityGradientDirectionReduced; - direction.GetBlock(gravity_displacement_force_test_utils::gravityPotentialValue) = 0.29; + direction.GetBlock(gravity_displacement_force_test_utils::gravityPotentialValue) = 0.29; - direction.GetBlock(gravity_displacement_force_test_utils::enthalpyValue) = -0.37; + direction.GetBlock(gravity_displacement_force_test_utils::enthalpyValue) = -0.37; direction.GetBlock(gravity_displacement_force_test_utils::barotropicConstantValue) = 0.43; @@ -673,7 +701,8 @@ TEST_CASE( mfem::Vector expectedDisplacementAction; preparedOperator.ApplyCompleteJacobianAction( - densityDirection, displacementDirection, gravityGradientDirection, expectedDisplacementAction + densityDirectionReduced, displacementDirectionReduced, gravityGradientDirectionReduced, + expectedDisplacementAction ); const mfem::Vector actualDisplacementAction = gravity_displacement_force_test_utils::copy_residual_block( @@ -707,4 +736,4 @@ TEST_CASE( for (const mfem::Vector &row : zeroRows) { CHECK(gravity_prepared_test_utils::global_norm(row, f.mesh->GetComm()) == 0.0); } -} \ No newline at end of file +} diff --git a/tests/operators/gravity_field.cpp b/tests/operators/gravity_field.cpp index cb7b087..be2ae4f 100644 --- a/tests/operators/gravity_field.cpp +++ b/tests/operators/gravity_field.cpp @@ -25,13 +25,18 @@ namespace { blocks::get_residual_block(blocks::gravity_field.poisson_term); blocks::form_layout make_gravity_layout(const fem::FEM &f) { + using DomainSchema = utils::domain::CoreEnvelopeVacuumDomainSchema; + const auto density_map = field::make_field_dof_map(*f.densityFes); + const auto displacement_map = field::make_field_dof_map(*f.displacementFes); + const auto flux_map = field::make_field_dof_map(*f.gravityFluxFes); + const auto potential_map = field::make_field_dof_map(*f.gravityPotentialFes); const std::array value_sizes{ - f.densityFes->GetTrueVSize(), f.displacementFes->GetTrueVSize(), f.gravityFluxFes->GetTrueVSize(), - f.gravityPotentialFes->GetTrueVSize() + density_map.reduced_size(), displacement_map.reduced_size(), flux_map.reduced_size(), + potential_map.reduced_size() }; const std::array residual_sizes{ - f.gravityFluxFes->GetTrueVSize(), f.gravityPotentialFes->GetTrueVSize() + flux_map.reduced_size(), potential_map.reduced_size() }; return blocks::form_layout(value_sizes, residual_sizes); @@ -623,16 +628,7 @@ namespace { using gravity_form = blocks::gravity_field_form; using gravity_layout = blocks::form_layout; gravity_layout make_gravity_jacobian_layout(const fem::FEM &f) { - const std::array value_sizes{ - f.densityFes->GetTrueVSize(), f.displacementFes->GetTrueVSize(), f.gravityFluxFes->GetTrueVSize(), - f.gravityPotentialFes->GetTrueVSize() - }; - - const std::array residual_sizes{ - f.gravityFluxFes->GetTrueVSize(), f.gravityPotentialFes->GetTrueVSize() - }; - - return gravity_layout(value_sizes, residual_sizes); + return make_gravity_layout(f); } template @@ -698,7 +694,9 @@ namespace { fill_value_block(state, layout, gravity_gradient_block, 0.4, 0.23); fill_value_block(state, layout, gravity_potential_block, 0.3, 0.37); - const mfem::Vector displacement = make_stateless_reference_displacement(f, deformed); + using DomainSchema = utils::domain::CoreEnvelopeVacuumDomainSchema; + const auto displacement_map = field::make_field_dof_map(*f.displacementFes); + const mfem::Vector displacement = displacement_map.gather(make_stateless_reference_displacement(f, deformed)); set_value_block(state, layout, displacement_block, displacement); return state; @@ -1026,15 +1024,23 @@ namespace { const mfem::Vector gravity_gradient = make_read_only_value_view(state, state_offsets, gravity_gradient_block); CHECK_THAT( - global_relative_vector_error(context.GetDensity(), density, communicator), + global_relative_vector_error( + context.GetDensityTrue(), context.GetDensityMap().scatter(density), communicator + ), Catch::Matchers::WithinAbs(0.0, 0.0) ); CHECK_THAT( - global_relative_vector_error(context.GetGeometryContext().GetDisplacement(), displacement, communicator), + global_relative_vector_error( + context.GetGeometryContext().GetDisplacementTrue(), context.GetDisplacementMap().scatter(displacement), + communicator + ), Catch::Matchers::WithinAbs(0.0, 0.0) ); CHECK_THAT( - global_relative_vector_error(context.GetGravityGradient(), gravity_gradient, communicator), + global_relative_vector_error( + context.GetGravityGradientTrue(), context.GetGravityGradientMap().scatter(gravity_gradient), + communicator + ), Catch::Matchers::WithinAbs(0.0, 0.0) ); } @@ -1179,7 +1185,7 @@ namespace { TEST_CASE( "Gravity Field Operator Mult Preserves Zero State And Static Blocks", - tags::unit &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_operator_unit ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -1256,7 +1262,7 @@ TEST_CASE( TEST_CASE( "Gravity Field Operator Mult Applies Stellar Source With Correct Sign", - tags::unit &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_operator_unit ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -1283,7 +1289,7 @@ TEST_CASE( .discretization = {1}, .displacement = {1}, .density = {1}, .gravity_gradient = {1}, .gravity_potential = {1} }; - const mfem::Vector density = make_constant_density(f, 1.0); + const mfem::Vector density = linearization_context.GetDensityMap().gather(make_constant_density(f, 1.0)); set_block(state, layout.value_offsets(), density_block, density); gravity_operator.Prepare(state, revisions); @@ -1319,7 +1325,7 @@ TEST_CASE( TEST_CASE( "Gravity Field Operator Mult Ignores Vacuum Density", - tags::unit &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_operator_unit ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -1347,8 +1353,9 @@ TEST_CASE( .discretization = {1}, .displacement = {1}, .density = {1}, .gravity_gradient = {1}, .gravity_potential = {1} }; - const mfem::Vector vacuum_density = make_vacuum_density(f, 1.0); - REQUIRE(vacuum_density.Norml2() > 0.0); + const mfem::Vector vacuum_density_true = make_vacuum_density(f, 1.0); + REQUIRE(vacuum_density_true.Norml2() > 0.0); + const mfem::Vector vacuum_density = linearization_context.GetDensityMap().gather(vacuum_density_true); set_block(state, layout.value_offsets(), density_block, vacuum_density); @@ -1362,7 +1369,7 @@ TEST_CASE( TEST_CASE( "Gravity Field Operator Mult Produces A Symmetric Positive Hdiv Mass " "Action", - tags::unit &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_operator_unit ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -1382,7 +1389,7 @@ TEST_CASE( f, *f.domainMapperStateless, linearization_context, layout.value_offsets(), gravity_jacobian ); - const mfem::Vector displacement = make_displacement(f); + const mfem::Vector displacement = linearization_context.GetDisplacementMap().gather(make_displacement(f)); const mfem::Vector gravity_gradient_a = make_test_vector(layout.size(gravity_gradient_block), 0.27); const mfem::Vector gravity_gradient_b = make_test_vector(layout.size(gravity_gradient_block), 1.13); @@ -1433,7 +1440,7 @@ TEST_CASE( TEST_CASE( "Gravity Field Operator Mult Is Additive In Fields At Fixed Displacement", - tags::unit &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_operator_unit ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -1453,8 +1460,8 @@ TEST_CASE( f, *f.domainMapperStateless, linearization_context, layout.value_offsets(), gravity_jacobian ); - const mfem::Vector displacement = make_displacement(f); - const mfem::Vector density = make_constant_density(f, 0.73); + const mfem::Vector displacement = linearization_context.GetDisplacementMap().gather(make_displacement(f)); + const mfem::Vector density = linearization_context.GetDensityMap().gather(make_constant_density(f, 0.73)); const mfem::Vector gravity_gradient = make_test_vector(layout.size(gravity_gradient_block), 0.41); const mfem::Vector gravity_potential = make_test_vector(layout.size(gravity_potential_block), 0.89); @@ -1514,7 +1521,7 @@ TEST_CASE( TEST_CASE( "Gravity Field Operator Stellar Hdiv Mass Matches Legacy Assembled Mass", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators &tags::legacy_comparison + tags::gravity_legacy ) { constexpr double parity_tolerance = 1.0e-10; constexpr double vacuum_tolerance = 1.0e-13; @@ -1631,7 +1638,7 @@ TEST_CASE( TEST_CASE( "Gravity Field Operator Hdiv Mass Action Matches Stateless Quadrature " "Reference", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_operator_integration ) { const utils::Args args = test_utils::setup_args(); @@ -1845,7 +1852,7 @@ TEST_CASE( TEST_CASE( "Gravity Field Jacobian Fixed Geometry Blocks Match Centered Differences", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_operator_integration ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -1865,9 +1872,10 @@ TEST_CASE( ); mfem::Vector state(layout.value_offsets().Last()); - state = 0.0; + state = 0.0; - const mfem::Vector displacement = make_stateless_reference_displacement(f, true); + const mfem::Vector displacement = + linearization_context.GetDisplacementMap().gather(make_stateless_reference_displacement(f, true)); set_value_block(state, layout, displacement_block, displacement); const operators::context::gravity_field::GravityFieldRevisions revisions{ @@ -1940,7 +1948,7 @@ TEST_CASE( TEST_CASE( "Gravity Field Jacobian Combined Fixed Geometry Direction Matches Centered " "Differences", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_operator_integration ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -1991,7 +1999,7 @@ TEST_CASE( TEST_CASE( "Gravity Field Jacobian Action Is Linear At Fixed Geometry", - tags::unit &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_operator_unit ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -2048,7 +2056,7 @@ TEST_CASE( } TEST_CASE( "Gravity Field Prepare Updates Shared Linearization Context", - tags::unit &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_operator_unit ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -2092,9 +2100,9 @@ TEST_CASE( CHECK(linearization_context.IsPrepared()); const double stored_state_norm = - global_vector_norm(linearization_context.GetDensity(), communicator) + - global_vector_norm(linearization_context.GetGeometryContext().GetDisplacement(), communicator) + - global_vector_norm(linearization_context.GetGravityGradient(), communicator); + global_vector_norm(linearization_context.GetDensityTrue(), communicator) + + global_vector_norm(linearization_context.GetGeometryContext().GetDisplacementTrue(), communicator) + + global_vector_norm(linearization_context.GetGravityGradientTrue(), communicator); CHECK(stored_state_norm > 0.0); @@ -2107,7 +2115,7 @@ TEST_CASE( } TEST_CASE( "Mapped Hdiv Mass Variation Is Linear In Displacement Direction", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_kernel_integration ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -2176,7 +2184,7 @@ TEST_CASE( TEST_CASE( "Mapped Hdiv Mass Variation Matches Centered Geometry Differences", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_kernel_integration ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -2220,7 +2228,7 @@ TEST_CASE( TEST_CASE( "Mapped Hdiv Mass Geometry Difference Converges To Analytic Variation", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators &tags::convergence + tags::gravity_kernel_convergence ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -2266,7 +2274,7 @@ TEST_CASE( TEST_CASE( "Mapped Source Variation Is Linear In Displacement Direction", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_kernel_integration ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -2339,7 +2347,7 @@ TEST_CASE( TEST_CASE( "Mapped Source Variation Matches Centered Geometry Differences", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_kernel_integration ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -2385,7 +2393,7 @@ TEST_CASE( TEST_CASE( "Mapped Source Geometry Difference Converges To Analytic Variation", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators &tags::convergence + tags::gravity_kernel_convergence ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -2431,7 +2439,7 @@ TEST_CASE( TEST_CASE( "Mapped Hdiv Mass Variation Is Symmetric", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_kernel_integration ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -2479,7 +2487,7 @@ TEST_CASE( TEST_CASE( "Mapped Geometry Variations Are Linear In Base Fields", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_kernel_integration ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -2571,7 +2579,7 @@ TEST_CASE( TEST_CASE( "Mapped Source Variation Ignores Vacuum Density", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_kernel_integration ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -2604,7 +2612,7 @@ TEST_CASE( TEST_CASE( "Gravity Field Jacobian Displacement Blocks Match Geometry Kernels", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_operator_integration ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -2618,9 +2626,6 @@ TEST_CASE( const HdivMassVariationTestFields fields = make_hdiv_mass_variation_test_fields(f); const mfem::Vector state = make_gravity_jacobian_state(f, layout, true); const mfem::Vector direction = make_displacement_direction(layout, fields.displacement_direction_1); - const mfem::Vector density = get_value_block(state, layout, density_block); - const mfem::Vector displacement = get_value_block(state, layout, displacement_block); - const mfem::Vector gravity_gradient = get_value_block(state, layout, gravity_gradient_block); operators::context::gravity_field::GravityFieldLinearizationContext linearization_context( f, *f.domainMapperStateless @@ -2645,19 +2650,30 @@ TEST_CASE( REQUIRE(jacobian_action.Size() == layout.residual_offsets().Last()); + const mfem::Vector &density_true = linearization_context.GetDensityTrue(); + const mfem::Vector &displacement_true = + linearization_context.GetGeometryContext().GetDisplacementTrue(); + const mfem::Vector &gravity_gradient_true = linearization_context.GetGravityGradientTrue(); + const mfem::Vector gradient_action = get_residual_block(jacobian_action, layout, gravity_gradient_residual_block); const mfem::Vector poisson_action = get_residual_block(jacobian_action, layout, gravity_poisson_residual_block); - mfem::Vector expected_gradient_action; - mfem::Vector expected_poisson_action; + mfem::Vector expected_gradient_action_true; + mfem::Vector expected_poisson_action_true; operators::kernels::apply_mapped_hdiv_mass_variation( - f, *f.domainMapperStateless, gravity_gradient, displacement, fields.displacement_direction_1, - expected_gradient_action + f, *f.domainMapperStateless, gravity_gradient_true, displacement_true, fields.displacement_direction_1, + expected_gradient_action_true ); operators::kernels::apply_mapped_source_variation( - f, *f.domainMapperStateless, density, displacement, fields.displacement_direction_1, expected_poisson_action + f, *f.domainMapperStateless, density_true, displacement_true, fields.displacement_direction_1, + expected_poisson_action_true ); - expected_poisson_action *= -1.0; + expected_poisson_action_true *= -1.0; + + const mfem::Vector expected_gradient_action = + linearization_context.GetGravityGradientMap().gather(expected_gradient_action_true); + const mfem::Vector expected_poisson_action = + linearization_context.GetGravityPotentialMap().gather(expected_poisson_action_true); MPI_Comm communicator = f.mesh->GetComm(); const double gradient_error = global_relative_error(gradient_action, expected_gradient_action, communicator); @@ -2677,7 +2693,7 @@ TEST_CASE( TEST_CASE( "Gravity Field Jacobian Displacement Direction Matches Centered " "Differences", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_operator_integration ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -2732,9 +2748,9 @@ TEST_CASE( linearization_context, state, layout.value_offsets(), communicator ); - const double base_density_norm = global_vector_norm(linearization_context.GetDensity(), communicator); + const double base_density_norm = global_vector_norm(linearization_context.GetDensityTrue(), communicator); const double base_gravity_gradient_norm = - global_vector_norm(linearization_context.GetGravityGradient(), communicator); + global_vector_norm(linearization_context.GetGravityGradientTrue(), communicator); INFO("Base density norm = " << base_density_norm); INFO("Base gravity-gradient norm = " << base_gravity_gradient_norm); @@ -2838,7 +2854,7 @@ TEST_CASE( } TEST_CASE( "Gravity Field Jacobian Displacement Difference Converges At Second Order", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators &tags::convergence + tags::gravity_operator_convergence ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -2964,7 +2980,7 @@ TEST_CASE( TEST_CASE( "Gravity Field Jacobian Complete Coupled Direction Matches Centered " "Differences", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_operator_integration ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -3108,7 +3124,7 @@ TEST_CASE( TEST_CASE( "Reduced Gravity Field Operator Solves Deformed Gravity System With MINRES", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators + tags::gravity_operator_integration ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -3118,11 +3134,11 @@ TEST_CASE( const gravity_layout layout = make_gravity_jacobian_layout(f); const HdivMassVariationTestFields fields = make_hdiv_mass_variation_test_fields(f); - const mfem::Vector density = make_source_variation_density(f); - const mfem::Vector displacement = fields.displacement; + const mfem::Vector density_true = make_source_variation_density(f); + const mfem::Vector displacement_true = fields.displacement; mfem::ParGridFunction legacy_displacement(f.displacementFes.get()); - legacy_displacement.SetFromTrueDofs(displacement); + legacy_displacement.SetFromTrueDofs(displacement_true); f.mapping->SetDisplacement(legacy_displacement); physics::update_stiffness_matrix(f); @@ -3133,6 +3149,8 @@ TEST_CASE( operators::context::gravity_field::GravityFieldLinearizationContext linearization_context( f, *f.domainMapperStateless ); + const mfem::Vector density = linearization_context.GetDensityMap().gather(density_true); + const mfem::Vector displacement = linearization_context.GetDisplacementMap().gather(displacement_true); operators::GravityFieldJacobianOperator gravity_jacobian( f, *f.domainMapperStateless, linearization_context, layout.value_offsets(), layout.residual_offsets() ); @@ -3291,7 +3309,7 @@ TEST_CASE( TEST_CASE( "Gravity Hdiv Operator Difference Is Localized To Vacuum Compactification", - tags::integration &tags::solver &tags::gravity &tags::mfem_operators &tags::legacy_comparison + tags::gravity_legacy ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -3311,16 +3329,7 @@ TEST_CASE( constexpr auto gravity_gradient_residual_block = mean_field::utils::blocks::get_residual_block(blocks::gravity_field.gradient_term); - const std::array value_sizes{ - f.densityFes->GetTrueVSize(), f.displacementFes->GetTrueVSize(), f.gravityFluxFes->GetTrueVSize(), - f.gravityPotentialFes->GetTrueVSize() - }; - - const std::array residual_sizes{ - f.gravityFluxFes->GetTrueVSize(), f.gravityPotentialFes->GetTrueVSize() - }; - - const blocks::form_layout layout(value_sizes, residual_sizes); + const blocks::form_layout layout = make_gravity_layout(f); const mfem::Vector gravity_gradient_true = make_full_support_gravity_gradient(f); diff --git a/tests/operators/prepared_displacement_operator.cpp b/tests/operators/prepared_displacement_operator.cpp index 3ddccb3..18dd46d 100644 --- a/tests/operators/prepared_displacement_operator.cpp +++ b/tests/operators/prepared_displacement_operator.cpp @@ -66,26 +66,24 @@ namespace prepared_displacement_residual_test_utils { } [[nodiscard]] mean_field::operators::DisplacementResidualLayout make_layout(const mean_field::fem::FEM &f) { + using DomainSchema = gravity_prepared_test_utils::DomainSchema; + + const auto densityMap = gravity_prepared_test_utils::make_field_map(f); + const auto displacementMap = gravity_prepared_test_utils::make_field_map(f); + const auto gravityFluxMap = + mean_field::field::make_field_dof_map(*f.gravityFluxFes); + const auto gravityPotentialMap = + mean_field::field::make_field_dof_map(*f.gravityPotentialFes); const auto enthalpyMap = make_enthalpy_map(f); - /* - * Transitional displacement-composer layout. - * - * Pressure has now migrated its h column to supported FieldDof - * coordinates, while the density-consuming mechanical children are - * intentionally still full-space until the next migration slice. - * - * The adapter only writes R_d. The unrelated residual-row sizes remain - * at their current full-space values in this standalone adapter test. - */ const std::array valueSizes{ - f.densityFes->GetTrueVSize(), f.displacementFes->GetTrueVSize(), f.gravityFluxFes->GetTrueVSize(), - f.gravityPotentialFes->GetTrueVSize(), enthalpyMap.reduced_size(), 1 + densityMap.reduced_size(), displacementMap.reduced_size(), gravityFluxMap.reduced_size(), + gravityPotentialMap.reduced_size(), enthalpyMap.reduced_size(), 1 }; const std::array residualSizes{ - f.gravityFluxFes->GetTrueVSize(), f.gravityPotentialFes->GetTrueVSize(), f.densityFes->GetTrueVSize(), - f.displacementFes->GetTrueVSize(), f.enthalpyFes->GetTrueVSize(), 1 + gravityFluxMap.reduced_size(), gravityPotentialMap.reduced_size(), densityMap.reduced_size(), + displacementMap.reduced_size(), enthalpyMap.reduced_size(), 1 }; return {valueSizes, residualSizes}; @@ -273,10 +271,10 @@ namespace prepared_displacement_residual_test_utils { const std::uint64_t potentialRevision ) { context.Prepare( - {.density = density, - .displacement = displacement, - .gravity_gradient = gravityGradient, - .gravity_potential = gravityPotential}, + {.density = context.GetDensityMap().gather(density), + .displacement = context.GetDisplacementMap().gather(displacement), + .gravity_gradient = context.GetGravityGradientMap().gather(gravityGradient), + .gravity_potential = context.GetGravityPotentialMap().gather(gravityPotential)}, make_gravity_revisions(dependencies, potentialRevision) ); } @@ -334,7 +332,7 @@ namespace prepared_displacement_residual_test_utils { TEST_CASE( "Prepared Displacement Residual Equals The Three Prepared Contributors", - tags::barotrope &tags::prepared &tags::integration &tags::residuals + tags::barotrope_prepared ) { using Operator = mean_field::operators::PreparedDisplacementResidualOperator; @@ -419,7 +417,7 @@ TEST_CASE( TEST_CASE( "Prepared Displacement Residual Selectively Orchestrates Its Children", - tags::barotrope &tags::prepared &tags::contexts &tags::integration + tags::barotrope_context_integration ) { mean_field::utils::Args args = test_utils::setup_args(); @@ -537,7 +535,7 @@ TEST_CASE( TEST_CASE( "Prepared Displacement Residual Jacobian Equals The Contributor Sums", - tags::barotrope &tags::prepared &tags::jacobian &tags::accuracy + tags::barotrope_prepared_jacobian_accuracy ) { mean_field::utils::Args args = test_utils::setup_args(); @@ -585,30 +583,37 @@ TEST_CASE( preparedOperator.Prepare({.enthalpy = enthalpy}, dependencies, rotation); + const mfem::Vector densityDirectionReduced = gravityContext.GetDensityMap().gather(densityDirection); + const mfem::Vector displacementDirectionReduced = gravityContext.GetDisplacementMap().gather(displacementDirection); + const mfem::Vector gravityDirectionReduced = gravityContext.GetGravityGradientMap().gather(gravityDirection); + mfem::Vector densityAction; mfem::Vector displacementAction; mfem::Vector gravityAction; mfem::Vector enthalpyAction; mfem::Vector completeAction; - preparedOperator.ApplyDensityJacobianAction(densityDirection, densityAction); + preparedOperator.ApplyDensityJacobianAction(densityDirectionReduced, densityAction); - preparedOperator.ApplyDisplacementJacobianAction(displacementDirection, displacementAction); + preparedOperator.ApplyDisplacementJacobianAction(displacementDirectionReduced, displacementAction); - preparedOperator.ApplyGravityGradientJacobianAction(gravityDirection, gravityAction); + preparedOperator.ApplyGravityGradientJacobianAction(gravityDirectionReduced, gravityAction); preparedOperator.ApplyEnthalpyJacobianAction(enthalpyDirection, enthalpyAction); preparedOperator.ApplyCompleteJacobianAction( - densityDirection, displacementDirection, gravityDirection, enthalpyDirection, completeAction + densityDirectionReduced, displacementDirectionReduced, gravityDirectionReduced, enthalpyDirection, + completeAction ); mfem::Vector expectedDensity; mfem::Vector expectedRotationDensity; - preparedOperator.GetGravityOperator().ApplyDensityJacobianAction(densityDirection, expectedDensity); + preparedOperator.GetGravityOperator().ApplyDensityJacobianAction(densityDirectionReduced, expectedDensity); - preparedOperator.GetRotationalOperator().ApplyDensityJacobianAction(densityDirection, expectedRotationDensity); + preparedOperator.GetRotationalOperator().ApplyDensityJacobianAction( + densityDirectionReduced, expectedRotationDensity + ); expectedDensity += expectedRotationDensity; @@ -616,21 +621,23 @@ TEST_CASE( mfem::Vector expectedGravityDisplacement; mfem::Vector expectedRotationDisplacement; - preparedOperator.GetPressureOperator().ApplyDisplacementJacobianAction(displacementDirection, expectedDisplacement); + preparedOperator.GetPressureOperator().ApplyDisplacementJacobianAction( + displacementDirectionReduced, expectedDisplacement + ); preparedOperator.GetGravityOperator().ApplyDisplacementJacobianAction( - displacementDirection, expectedGravityDisplacement + displacementDirectionReduced, expectedGravityDisplacement ); preparedOperator.GetRotationalOperator().ApplyDisplacementJacobianAction( - displacementDirection, expectedRotationDisplacement + displacementDirectionReduced, expectedRotationDisplacement ); expectedDisplacement += expectedGravityDisplacement; expectedDisplacement += expectedRotationDisplacement; mfem::Vector expectedGravity; - preparedOperator.GetGravityOperator().ApplyGravityGradientJacobianAction(gravityDirection, expectedGravity); + preparedOperator.GetGravityOperator().ApplyGravityGradientJacobianAction(gravityDirectionReduced, expectedGravity); mfem::Vector expectedEnthalpy; preparedOperator.GetPressureOperator().ApplyEnthalpyJacobianAction(enthalpyDirection, expectedEnthalpy); @@ -672,7 +679,7 @@ TEST_CASE( TEST_CASE( "Prepared Displacement Residual Jacobian Matches A Simultaneous " "Centered Difference On Deformed Geometry", - tags::barotrope &tags::prepared &tags::jacobian &tags::accuracy &tags::geometry + tags::barotrope_prepared_jacobian_accuracy &tags::geometry ) { mean_field::utils::Args args = test_utils::setup_args(); @@ -720,10 +727,15 @@ TEST_CASE( preparedOperator.Prepare({.enthalpy = baseEnthalpy}, dependencies, rotation); + const mfem::Vector densityDirectionReduced = gravityContext.GetDensityMap().gather(densityDirection); + const mfem::Vector displacementDirectionReduced = gravityContext.GetDisplacementMap().gather(displacementDirection); + const mfem::Vector gravityDirectionReduced = gravityContext.GetGravityGradientMap().gather(gravityDirection); + mfem::Vector jacobianAction; preparedOperator.ApplyCompleteJacobianAction( - densityDirection, displacementDirection, gravityDirection, enthalpyDirection, jacobianAction + densityDirectionReduced, displacementDirectionReduced, gravityDirectionReduced, enthalpyDirection, + jacobianAction ); constexpr double step = 1.0e-5; @@ -793,7 +805,7 @@ TEST_CASE( TEST_CASE( "Prepared Displacement Residual MFEM Adapter Routes Only R-d", - tags::barotrope &tags::prepared &tags::jacobian &tags::mfem_operators &tags::unit + tags::barotrope_prepared_jacobian_unit ) { using JacobianForm = mean_field::utils::blocks::barotropic_equilibrium_jacobian_form; @@ -876,6 +888,10 @@ TEST_CASE( preparedOperator.Prepare({.enthalpy = enthalpy}, dependencies, rotation); + const mfem::Vector densityDirectionReduced = gravityContext.GetDensityMap().gather(densityDirection); + const mfem::Vector displacementDirectionReduced = gravityContext.GetDisplacementMap().gather(displacementDirection); + const mfem::Vector gravityDirectionReduced = gravityContext.GetGravityGradientMap().gather(gravityDirection); + const mean_field::operators::DisplacementResidualLayout layout = prepared_displacement_residual_test_utils::make_layout(f); @@ -889,24 +905,25 @@ TEST_CASE( ); mfem::BlockVector direction(layout.value_offsets()); - direction = 0.0; + direction = 0.0; - direction.GetBlock(prepared_displacement_residual_test_utils::densityValue) = densityDirection; + direction.GetBlock(prepared_displacement_residual_test_utils::densityValue) = densityDirectionReduced; - direction.GetBlock(prepared_displacement_residual_test_utils::displacementValue) = displacementDirection; + direction.GetBlock(prepared_displacement_residual_test_utils::displacementValue) = displacementDirectionReduced; - direction.GetBlock(prepared_displacement_residual_test_utils::gravityGradientValue) = gravityDirection; + direction.GetBlock(prepared_displacement_residual_test_utils::gravityGradientValue) = gravityDirectionReduced; - direction.GetBlock(prepared_displacement_residual_test_utils::enthalpyValue) = enthalpyDirection; + direction.GetBlock(prepared_displacement_residual_test_utils::enthalpyValue) = enthalpyDirection; - direction.GetBlock(prepared_displacement_residual_test_utils::gravityPotentialValue) = 0.59; + direction.GetBlock(prepared_displacement_residual_test_utils::gravityPotentialValue) = 0.59; direction.GetBlock(prepared_displacement_residual_test_utils::barotropicConstantValue) = -0.73; mfem::Vector expectedDisplacementAction; preparedOperator.ApplyCompleteJacobianAction( - densityDirection, displacementDirection, gravityDirection, enthalpyDirection, expectedDisplacementAction + densityDirectionReduced, displacementDirectionReduced, gravityDirectionReduced, enthalpyDirection, + expectedDisplacementAction ); mfem::Vector action; diff --git a/tests/operators/prepared_gravity_source.cpp b/tests/operators/prepared_gravity_source.cpp index 2060d2f..84a0f10 100644 --- a/tests/operators/prepared_gravity_source.cpp +++ b/tests/operators/prepared_gravity_source.cpp @@ -11,31 +11,36 @@ namespace prepared_test = gravity_prepared_test_utils; TEST_CASE( "Prepared Mapped Gravity Source Matches Stateless Kernel", - tags::integration &tags::mfem_operators &tags::prepared + tags::gravity_prepared ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); operators::PreparedMappedGravitySourceOperator prepared_operator(f, *f.domainMapperStateless); - REQUIRE(prepared_operator.Width() == f.densityFes->GetTrueVSize()); - REQUIRE(prepared_operator.Height() == f.gravityPotentialFes->GetTrueVSize()); + REQUIRE(prepared_operator.Width() == prepared_operator.GetDensityMap().reduced_size()); + REQUIRE(prepared_operator.Height() == prepared_operator.GetPotentialMap().reduced_size()); - const mfem::Vector density = prepared_test::make_deterministic_vector(f.densityFes->GetTrueVSize(), 0.41); - const MPI_Comm communicator = f.mesh->GetComm(); + const mfem::Vector density_true = prepared_test::make_deterministic_vector(f.densityFes->GetTrueVSize(), 0.41); + const mfem::Vector density = prepared_operator.GetDensityMap().gather(density_true); + const MPI_Comm communicator = f.mesh->GetComm(); mfem::Vector identity_action; mfem::Vector deformed_action; for (const double deformation_scale : {0.0, 1.0}) { - const mfem::Vector displacement = prepared_test::make_displacement(f, deformation_scale); + const mfem::Vector displacement_true = prepared_test::make_displacement(f, deformation_scale); + const mfem::Vector displacement = prepared_operator.GetDisplacementMap().gather(displacement_true); prepared_operator.Prepare(displacement); mfem::Vector prepared_action; - mfem::Vector reference_action; prepared_operator.Mult(density, prepared_action); - operators::kernels::apply_mapped_source(f, *f.domainMapperStateless, density, displacement, reference_action); + mfem::Vector reference_action_true; + operators::kernels::apply_mapped_source( + f, *f.domainMapperStateless, density_true, displacement_true, reference_action_true + ); + const mfem::Vector reference_action = prepared_operator.GetPotentialMap().gather(reference_action_true); const double relative_error = prepared_test::relative_error(prepared_action, reference_action, communicator); @@ -64,28 +69,36 @@ TEST_CASE( TEST_CASE( "Prepared Mapped Gravity Source Preserves Linearity And Excludes Vacuum", - tags::integration &tags::gravity &tags::prepared + tags::gravity_prepared ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); operators::PreparedMappedGravitySourceOperator prepared_operator(f, *f.domainMapperStateless); - REQUIRE(prepared_operator.Width() == f.densityFes->GetTrueVSize()); - REQUIRE(prepared_operator.Height() == f.gravityPotentialFes->GetTrueVSize()); - const mfem::Vector displacement = prepared_test::make_displacement(f, 1.0); + REQUIRE(prepared_operator.Width() == prepared_operator.GetDensityMap().reduced_size()); + REQUIRE(prepared_operator.Height() == prepared_operator.GetPotentialMap().reduced_size()); + const mfem::Vector displacement = + prepared_operator.GetDisplacementMap().gather(prepared_test::make_displacement(f, 1.0)); prepared_operator.Prepare(displacement); - const mfem::Vector first = prepared_test::make_deterministic_vector(f.densityFes->GetTrueVSize(), 0.27); - const mfem::Vector second = prepared_test::make_deterministic_vector(f.densityFes->GetTrueVSize(), 0.79); - const mfem::Vector combination = prepared_test::linear_combination(first, 1.3, second, -0.6); - const mfem::Vector stellar_density = prepared_test::make_domain_supported_density(f, true); - const mfem::Vector vacuum_density = prepared_test::make_domain_supported_density(f, false); + const mfem::Vector first = prepared_operator.GetDensityMap().gather( + prepared_test::make_deterministic_vector(f.densityFes->GetTrueVSize(), 0.27) + ); + const mfem::Vector second = prepared_operator.GetDensityMap().gather( + prepared_test::make_deterministic_vector(f.densityFes->GetTrueVSize(), 0.79) + ); + const mfem::Vector combination = prepared_test::linear_combination(first, 1.3, second, -0.6); + const mfem::Vector stellar_density = + prepared_operator.GetDensityMap().gather(prepared_test::make_domain_supported_density(f, true)); + const mfem::Vector vacuum_density = + prepared_operator.GetDensityMap().gather(prepared_test::make_domain_supported_density(f, false)); mfem::Vector first_action; mfem::Vector second_action; mfem::Vector combination_action; mfem::Vector stellar_action; mfem::Vector vacuum_action; + mfem::Vector transpose_action; prepared_operator.Mult(first, first_action); prepared_operator.Mult(second, second_action); @@ -93,6 +106,10 @@ TEST_CASE( prepared_operator.Mult(stellar_density, stellar_action); prepared_operator.Mult(vacuum_density, vacuum_action); + const mfem::Vector potential = + prepared_test::make_deterministic_vector(prepared_operator.GetPotentialMap().reduced_size(), 1.17); + prepared_operator.MultTranspose(potential, transpose_action); + const mfem::Vector expected_combination = prepared_test::linear_combination(first_action, 1.3, second_action, -0.6); const MPI_Comm communicator = f.mesh->GetComm(); @@ -100,6 +117,9 @@ TEST_CASE( prepared_test::relative_error(combination_action, expected_combination, communicator); const double stellar_norm = prepared_test::global_norm(stellar_action, communicator); const double vacuum_norm = prepared_test::global_norm(vacuum_action, communicator); + const double forward_pairing = prepared_test::global_dot(first_action, potential, communicator); + const double transpose_pairing = prepared_test::global_dot(first, transpose_action, communicator); + const double transpose_error = prepared_test::relative_scalar_error(forward_pairing, transpose_pairing); const std::uint64_t preparation_count = prepared_operator.GetPreparationCount(); mfem::Vector repeated_action; @@ -108,10 +128,12 @@ TEST_CASE( INFO("Relative source linearity error = " << linearity_error); INFO("Stellar source norm = " << stellar_norm); INFO("Vacuum-only source norm = " << vacuum_norm); + INFO("Relative source-transpose pairing error = " << transpose_error); CHECK_THAT(linearity_error, WithinAbs(0.0, 2.0e-12)); CHECK(stellar_norm > 0.0); CHECK(vacuum_norm <= 1.0e-13 * stellar_norm); + CHECK_THAT(transpose_error, WithinAbs(0.0, 2.0e-12)); CHECK(prepared_test::relative_error(repeated_action, first_action, communicator) < 2.0e-14); CHECK(prepared_operator.GetPreparationCount() == preparation_count); } diff --git a/tests/operators/prepared_hdiv_mass.cpp b/tests/operators/prepared_hdiv_mass.cpp index 57e54ad..ee8b586 100644 --- a/tests/operators/prepared_hdiv_mass.cpp +++ b/tests/operators/prepared_hdiv_mass.cpp @@ -11,32 +11,37 @@ namespace prepared_test = gravity_prepared_test_utils; TEST_CASE( "Prepared Mapped Hdiv Mass Matches Stateless Kernel", - tags::integration &tags::mfem_operators &tags::prepared + tags::gravity_prepared ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); operators::PreparedMappedHDivMassOperator prepared_operator(f, *f.domainMapperStateless); + REQUIRE(prepared_operator.Width() == prepared_operator.GetFluxMap().reduced_size()); + REQUIRE(prepared_operator.Height() == prepared_operator.GetFluxMap().reduced_size()); - const mfem::Vector gravity_gradient = + const mfem::Vector gravity_gradient_true = prepared_test::make_deterministic_vector(f.gravityFluxFes->GetTrueVSize(), 0.21); - const MPI_Comm communicator = f.gravityFluxFes->GetComm(); + const mfem::Vector gravity_gradient = prepared_operator.GetFluxMap().gather(gravity_gradient_true); + const MPI_Comm communicator = f.gravityFluxFes->GetComm(); mfem::Vector identity_action; mfem::Vector deformed_action; for (const double deformation_scale : {0.0, 1.0}) { - const mfem::Vector displacement = prepared_test::make_displacement(f, deformation_scale); + const mfem::Vector displacement_true = prepared_test::make_displacement(f, deformation_scale); + const mfem::Vector displacement = prepared_operator.GetDisplacementMap().gather(displacement_true); prepared_operator.Prepare(displacement); mfem::Vector prepared_action; - mfem::Vector reference_action; prepared_operator.Mult(gravity_gradient, prepared_action); + mfem::Vector reference_action_true; operators::kernels::apply_mapped_hdiv_mass( - f, *f.domainMapperStateless, gravity_gradient, displacement, reference_action + f, *f.domainMapperStateless, gravity_gradient_true, displacement_true, reference_action_true ); + const mfem::Vector reference_action = prepared_operator.GetFluxMap().gather(reference_action_true); const double relative_error = prepared_test::relative_error(prepared_action, reference_action, communicator); @@ -65,17 +70,24 @@ TEST_CASE( TEST_CASE( "Prepared Mapped Hdiv Mass Preserves Operator Identities", - tags::integration &tags::gravity &tags::prepared + tags::gravity_prepared ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); operators::PreparedMappedHDivMassOperator prepared_operator(f, *f.domainMapperStateless); - const mfem::Vector displacement = prepared_test::make_displacement(f, 1.0); + REQUIRE(prepared_operator.Width() == prepared_operator.GetFluxMap().reduced_size()); + REQUIRE(prepared_operator.Height() == prepared_operator.GetFluxMap().reduced_size()); + const mfem::Vector displacement = + prepared_operator.GetDisplacementMap().gather(prepared_test::make_displacement(f, 1.0)); prepared_operator.Prepare(displacement); - const mfem::Vector first = prepared_test::make_deterministic_vector(f.gravityFluxFes->GetTrueVSize(), 0.17); - const mfem::Vector second = prepared_test::make_deterministic_vector(f.gravityFluxFes->GetTrueVSize(), 0.83); + const mfem::Vector first = prepared_operator.GetFluxMap().gather( + prepared_test::make_deterministic_vector(f.gravityFluxFes->GetTrueVSize(), 0.17) + ); + const mfem::Vector second = prepared_operator.GetFluxMap().gather( + prepared_test::make_deterministic_vector(f.gravityFluxFes->GetTrueVSize(), 0.83) + ); const mfem::Vector combination = prepared_test::linear_combination(first, 1.7, second, -0.4); mfem::Vector first_action; @@ -121,4 +133,4 @@ TEST_CASE( CHECK(second_energy > 0.0); CHECK(prepared_test::relative_error(repeated_action, first_action, communicator) < 2.0e-14); CHECK(prepared_operator.GetPreparationCount() == preparation_count); -} \ No newline at end of file +} diff --git a/tests/operators/prepared_hydrostatic_equilibrium.cpp b/tests/operators/prepared_hydrostatic_equilibrium.cpp index f4b926e..0a96bab 100644 --- a/tests/operators/prepared_hydrostatic_equilibrium.cpp +++ b/tests/operators/prepared_hydrostatic_equilibrium.cpp @@ -34,14 +34,18 @@ namespace prepared_hydrostatic_test_utils { const mean_field::fem::FEM &f, const double phase = 0.19 ) { - return gravity_prepared_test_utils::make_deterministic_vector(f.enthalpyFes->GetTrueVSize(), phase); + return field_dof_test_utils::make_deterministic_supported_vector( + *f.enthalpyFes, phase + ); } mfem::Vector make_gravity_potential( const mean_field::fem::FEM &f, const double phase = 0.37 ) { - return gravity_prepared_test_utils::make_deterministic_vector(f.gravityPotentialFes->GetTrueVSize(), phase); + return field_dof_test_utils::make_deterministic_supported_vector( + *f.gravityPotentialFes, phase + ); } mean_field::physics::RigidRotation make_rotation(const double scale = 1.0) { @@ -63,7 +67,7 @@ namespace prepared_hydrostatic_test_utils { TEST_CASE( "Prepared Hydrostatic Residual Matches Stateless Kernel", - tags::barotrope &tags::hydro &tags::prepared &tags::residuals &tags::unit + tags::barotrope_hydrostatic_prepared_residual &tags::unit ) { auto args = test_utils::setup_args(); @@ -75,7 +79,7 @@ TEST_CASE( const mfem::Vector gravityPotential = prepared_hydrostatic_test_utils::make_gravity_potential(f); - const mfem::Vector displacement = gravity_prepared_test_utils::make_displacement(f, 0.73); + const mfem::Vector displacement = field_dof_test_utils::make_supported_displacement(f, 0.73); constexpr double bernoulliConstant = 0.41; @@ -97,9 +101,8 @@ TEST_CASE( preparedOperator.BuildResidual(preparedResidual); - mean_field::operators::kernels::apply_hydrostatic_equilibrium( - f, *f.domainMapperStateless, rotation, enthalpy, gravityPotential, displacement, bernoulliConstant, - referenceResidual + field_dof_test_utils::apply_hydrostatic_reference( + f, rotation, enthalpy, gravityPotential, displacement, bernoulliConstant, referenceResidual ); const double relativeError = @@ -131,7 +134,7 @@ TEST_CASE( TEST_CASE( "Prepared Hydrostatic Residual Reuses And Selectively Rebuilds Data", - tags::barotrope &tags::hydro &tags::prepared &tags::residuals &tags::unit + tags::barotrope_hydrostatic_prepared_residual &tags::unit ) { auto args = test_utils::setup_args(); @@ -143,7 +146,7 @@ TEST_CASE( mfem::Vector gravityPotential = prepared_hydrostatic_test_utils::make_gravity_potential(f); - mfem::Vector displacement = gravity_prepared_test_utils::make_displacement(f, 0.42); + mfem::Vector displacement = field_dof_test_utils::make_supported_displacement(f, 0.42); double bernoulliConstant = 0.36; @@ -168,7 +171,7 @@ TEST_CASE( enthalpy(0) += 0.29; gravityPotential(0) -= 0.17; - displacement = gravity_prepared_test_utils::make_displacement(f, 0.73); + displacement = field_dof_test_utils::make_supported_displacement(f, 0.73); bernoulliConstant += 0.23; const auto unchangedReport = preparedOperator.Prepare( @@ -199,9 +202,9 @@ TEST_CASE( preparedOperator.BuildResidual(enthalpyResidual); - mean_field::operators::kernels::apply_hydrostatic_equilibrium( - f, *f.domainMapperStateless, initialRotation, enthalpy, initialGravityPotential, initialDisplacement, - initialBernoulliConstant, enthalpyReference + field_dof_test_utils::apply_hydrostatic_reference( + f, initialRotation, enthalpy, initialGravityPotential, initialDisplacement, initialBernoulliConstant, + enthalpyReference ); CHECK_FALSE(enthalpyReport.contextReport.preparedStaticDependencies); @@ -226,9 +229,9 @@ TEST_CASE( preparedOperator.BuildResidual(rotationResidual); - mean_field::operators::kernels::apply_hydrostatic_equilibrium( - f, *f.domainMapperStateless, changedRotation, enthalpy, initialGravityPotential, initialDisplacement, - initialBernoulliConstant, rotationReference + field_dof_test_utils::apply_hydrostatic_reference( + f, changedRotation, enthalpy, initialGravityPotential, initialDisplacement, initialBernoulliConstant, + rotationReference ); CHECK_FALSE(rotationReport.contextReport.preparedStaticDependencies); @@ -253,9 +256,9 @@ TEST_CASE( preparedOperator.BuildResidual(displacementResidual); - mean_field::operators::kernels::apply_hydrostatic_equilibrium( - f, *f.domainMapperStateless, changedRotation, enthalpy, initialGravityPotential, displacement, - initialBernoulliConstant, displacementReference + field_dof_test_utils::apply_hydrostatic_reference( + f, changedRotation, enthalpy, initialGravityPotential, displacement, initialBernoulliConstant, + displacementReference ); CHECK_FALSE(displacementReport.contextReport.preparedStaticDependencies); diff --git a/tests/operators/prepared_hydrostatic_equilibrium_analytic_accuracy.cpp b/tests/operators/prepared_hydrostatic_equilibrium_analytic_accuracy.cpp index 73906ce..0f9d4a5 100644 --- a/tests/operators/prepared_hydrostatic_equilibrium_analytic_accuracy.cpp +++ b/tests/operators/prepared_hydrostatic_equilibrium_analytic_accuracy.cpp @@ -144,8 +144,7 @@ namespace prepared_hydrostatic_analytic_solve_test_utils { TEST_CASE( "Prepared Hydrostatic Operator Solves Analytic Bernoulli Equilibria", - tags::barotrope &tags::hydro &tags::prepared &tags::integration &tags::solver &tags::convergence &tags::accuracy - &tags::analytic_comparison + tags::barotrope_hydrostatic_prepared_analytic &tags::convergence &tags::accuracy ) { using prepared_hydrostatic_analytic_solve_test_utils::AnalyticCase; @@ -181,6 +180,15 @@ TEST_CASE( const MPI_Comm communicator = f.mesh->GetComm(); + const mean_field::field::FieldDofMap enthalpyMap = + field_dof_test_utils::make_map(*f.enthalpyFes); + + const mean_field::field::FieldDofMap gravityPotentialMap = + field_dof_test_utils::make_map(*f.gravityPotentialFes); + + const mean_field::field::FieldDofMap displacementMap = + field_dof_test_utils::make_map(*f.displacementFes); + const mfem::Array stellarElementMarker = prepared_hydrostatic_analytic_solve_test_utils::make_stellar_element_marker(f); @@ -234,11 +242,14 @@ TEST_CASE( potentialField.ProjectCoefficient(potentialCoefficient); - mfem::Vector displacement; - mfem::Vector gravityPotential; + mfem::Vector displacementTrue; + mfem::Vector gravityPotentialTrue; - displacementField.GetTrueDofs(displacement); - potentialField.GetTrueDofs(gravityPotential); + displacementField.GetTrueDofs(displacementTrue); + potentialField.GetTrueDofs(gravityPotentialTrue); + + const mfem::Vector displacement = displacementMap.gather(displacementTrue); + const mfem::Vector gravityPotential = gravityPotentialMap.gather(gravityPotentialTrue); /* * This projection is not used as the solution. It gives @@ -266,7 +277,7 @@ TEST_CASE( /* * Begin deliberately far from equilibrium. */ - mfem::Vector enthalpy(f.enthalpyFes->GetTrueVSize()); + mfem::Vector enthalpy(enthalpyMap.reduced_size()); enthalpy = 0.0; @@ -301,21 +312,20 @@ TEST_CASE( * enthalpy solve. */ prepared_hydrostatic_analytic_solve_test_utils::EnthalpyJacobianOperator enthalpyJacobian( - f.enthalpyFes->GetTrueVSize(), preparedOperator + enthalpyMap.reduced_size(), preparedOperator ); mfem::Vector rightHandSide(initialResidual); rightHandSide *= -1.0; - mfem::Vector enthalpyCorrection(f.enthalpyFes->GetTrueVSize()); + mfem::Vector enthalpyCorrection(enthalpyMap.reduced_size()); enthalpyCorrection = 0.0; /* - * The operator is positive definite on stellar-supported - * enthalpy DOFs and semidefinite on exterior-only DOFs. - * The RHS is in its range, so MINRES is appropriate for - * the compatible system. + * The reduced operator contains only stellar-supported + * enthalpy DOFs and is positive definite. MINRES remains + * appropriate for this symmetric system. */ mfem::MINRESSolver linearSolver(communicator); @@ -375,7 +385,9 @@ TEST_CASE( */ mfem::ParGridFunction solvedEnthalpyField(f.enthalpyFes.get()); - solvedEnthalpyField.SetFromTrueDofs(enthalpy); + mfem::Vector enthalpyTrue(enthalpyMap.full_size()); + enthalpyMap.scatter(enthalpy, enthalpyTrue); + solvedEnthalpyField.SetFromTrueDofs(enthalpyTrue); const double solvedAnalyticError = solvedEnthalpyField.ComputeL2Error(exactEnthalpyCoefficient, nullptr, &stellarElementMarker); @@ -416,4 +428,4 @@ TEST_CASE( CHECK(relativeSolvedAnalyticError / relativeProjectionError < 5.0); } } -} \ No newline at end of file +} diff --git a/tests/operators/prepared_hydrostatic_equilibrium_complete_jacobian.cpp b/tests/operators/prepared_hydrostatic_equilibrium_complete_jacobian.cpp index 8482ea9..34c9277 100644 --- a/tests/operators/prepared_hydrostatic_equilibrium_complete_jacobian.cpp +++ b/tests/operators/prepared_hydrostatic_equilibrium_complete_jacobian.cpp @@ -51,9 +51,9 @@ namespace prepared_hydrostatic_complete_test_utils { const double firstPhase, const double secondPhase ) { - mfem::Vector direction = gravity_prepared_test_utils::make_displacement(f, firstPhase); + mfem::Vector direction = field_dof_test_utils::make_supported_displacement(f, firstPhase); - const mfem::Vector secondField = gravity_prepared_test_utils::make_displacement(f, secondPhase); + const mfem::Vector secondField = field_dof_test_utils::make_supported_displacement(f, secondPhase); direction -= secondField; return direction; @@ -97,14 +97,12 @@ namespace prepared_hydrostatic_complete_test_utils { mfem::Vector residualPlus; mfem::Vector residualMinus; - mean_field::operators::kernels::apply_hydrostatic_equilibrium( - f, *f.domainMapperStateless, rotation, enthalpyPlus, gravityPotentialPlus, displacementPlus, - bernoulliConstantPlus, residualPlus + field_dof_test_utils::apply_hydrostatic_reference( + f, rotation, enthalpyPlus, gravityPotentialPlus, displacementPlus, bernoulliConstantPlus, residualPlus ); - mean_field::operators::kernels::apply_hydrostatic_equilibrium( - f, *f.domainMapperStateless, rotation, enthalpyMinus, gravityPotentialMinus, displacementMinus, - bernoulliConstantMinus, residualMinus + field_dof_test_utils::apply_hydrostatic_reference( + f, rotation, enthalpyMinus, gravityPotentialMinus, displacementMinus, bernoulliConstantMinus, residualMinus ); difference = residualPlus; @@ -131,7 +129,7 @@ namespace prepared_hydrostatic_complete_test_utils { TEST_CASE( "Prepared Hydrostatic Complete Jacobian Matches Sum And Centered " "Differences", - tags::barotrope &tags::hydro &tags::integration &tags::jacobian &tags::prepared &tags::self_consistency + tags::barotrope_hydrostatic_prepared_jacobian &tags::self_consistency ) { auto args = test_utils::setup_args(); @@ -140,12 +138,16 @@ TEST_CASE( mean_field::operators::PreparedHydrostaticEquilibriumOperator preparedOperator(f, *f.domainMapperStateless); const mfem::Vector enthalpy = - gravity_prepared_test_utils::make_deterministic_vector(f.enthalpyFes->GetTrueVSize(), 0.34); + field_dof_test_utils::make_deterministic_supported_vector( + *f.enthalpyFes, 0.34 + ); const mfem::Vector gravityPotential = - gravity_prepared_test_utils::make_deterministic_vector(f.gravityPotentialFes->GetTrueVSize(), 0.57); + field_dof_test_utils::make_deterministic_supported_vector( + *f.gravityPotentialFes, 0.57 + ); - const mfem::Vector displacement = gravity_prepared_test_utils::make_displacement(f, 0.68); + const mfem::Vector displacement = field_dof_test_utils::make_supported_displacement(f, 0.68); constexpr double bernoulliConstant = 0.43; @@ -159,10 +161,14 @@ TEST_CASE( ); const mfem::Vector enthalpyVariation = - gravity_prepared_test_utils::make_deterministic_vector(f.enthalpyFes->GetTrueVSize(), 1.07); + field_dof_test_utils::make_deterministic_supported_vector( + *f.enthalpyFes, 1.07 + ); const mfem::Vector gravityPotentialVariation = - gravity_prepared_test_utils::make_deterministic_vector(f.gravityPotentialFes->GetTrueVSize(), 1.31); + field_dof_test_utils::make_deterministic_supported_vector( + *f.gravityPotentialFes, 1.31 + ); const mfem::Vector displacementVariation = prepared_hydrostatic_complete_test_utils::make_displacement_direction(f, 1.19, 0.38); @@ -221,7 +227,7 @@ TEST_CASE( TEST_CASE( "Prepared Hydrostatic MFEM Adapter Uses Four Block Layout And Reuses " "Preparation", - tags::barotrope &tags::hydro &tags::integration &tags::jacobian &tags::mfem_operators &tags::prepared &tags::unit + tags::barotrope_hydrostatic_prepared_jacobian &tags::mfem_operators &tags::unit ) { auto args = test_utils::setup_args(); @@ -230,12 +236,16 @@ TEST_CASE( mean_field::operators::PreparedHydrostaticEquilibriumOperator preparedOperator(f, *f.domainMapperStateless); const mfem::Vector enthalpy = - gravity_prepared_test_utils::make_deterministic_vector(f.enthalpyFes->GetTrueVSize(), 0.41); + field_dof_test_utils::make_deterministic_supported_vector( + *f.enthalpyFes, 0.41 + ); const mfem::Vector gravityPotential = - gravity_prepared_test_utils::make_deterministic_vector(f.gravityPotentialFes->GetTrueVSize(), 0.63); + field_dof_test_utils::make_deterministic_supported_vector( + *f.gravityPotentialFes, 0.63 + ); - const mfem::Vector displacement = gravity_prepared_test_utils::make_displacement(f, 0.74); + const mfem::Vector displacement = field_dof_test_utils::make_supported_displacement(f, 0.74); constexpr double bernoulliConstant = 0.38; @@ -255,20 +265,26 @@ TEST_CASE( CHECK( layout.Offset(mean_field::operators::HydrostaticJacobianInputBlock::gravityPotential) == - f.enthalpyFes->GetTrueVSize() + preparedOperator.GetEnthalpyMap().reduced_size() ); CHECK(layout.Size(mean_field::operators::HydrostaticJacobianInputBlock::bernoulliConstant) == 1); CHECK(adapter.Width() == layout.GetTotalSize()); CHECK(adapter.Height() == layout.GetResidualSize()); - CHECK(adapter.Height() == f.enthalpyFes->GetTrueVSize()); + CHECK(adapter.Height() == preparedOperator.GetEnthalpyMap().reduced_size()); + + CHECK(preparedOperator.GetEnthalpyMap().inactive_size() > 0); const mfem::Vector enthalpyVariation = - gravity_prepared_test_utils::make_deterministic_vector(f.enthalpyFes->GetTrueVSize(), 1.12); + field_dof_test_utils::make_deterministic_supported_vector( + *f.enthalpyFes, 1.12 + ); const mfem::Vector gravityPotentialVariation = - gravity_prepared_test_utils::make_deterministic_vector(f.gravityPotentialFes->GetTrueVSize(), 1.39); + field_dof_test_utils::make_deterministic_supported_vector( + *f.gravityPotentialFes, 1.39 + ); const mfem::Vector displacementVariation = prepared_hydrostatic_complete_test_utils::make_displacement_direction(f, 1.28, 0.49); diff --git a/tests/operators/prepared_hydrostatic_equilibrium_displacement_jacobian.cpp b/tests/operators/prepared_hydrostatic_equilibrium_displacement_jacobian.cpp index 38c5a7e..4436ceb 100644 --- a/tests/operators/prepared_hydrostatic_equilibrium_displacement_jacobian.cpp +++ b/tests/operators/prepared_hydrostatic_equilibrium_displacement_jacobian.cpp @@ -34,14 +34,18 @@ namespace prepared_hydrostatic_displacement_test_utils { const mean_field::fem::FEM &f, const double phase = 0.29 ) { - return gravity_prepared_test_utils::make_deterministic_vector(f.enthalpyFes->GetTrueVSize(), phase); + return field_dof_test_utils::make_deterministic_supported_vector( + *f.enthalpyFes, phase + ); } mfem::Vector make_gravity_potential( const mean_field::fem::FEM &f, const double phase = 0.47 ) { - return gravity_prepared_test_utils::make_deterministic_vector(f.gravityPotentialFes->GetTrueVSize(), phase); + return field_dof_test_utils::make_deterministic_supported_vector( + *f.gravityPotentialFes, phase + ); } mfem::Vector make_displacement_direction( @@ -49,9 +53,9 @@ namespace prepared_hydrostatic_displacement_test_utils { const double firstPhase, const double secondPhase ) { - mfem::Vector direction = gravity_prepared_test_utils::make_displacement(f, firstPhase); + mfem::Vector direction = field_dof_test_utils::make_supported_displacement(f, firstPhase); - const mfem::Vector secondField = gravity_prepared_test_utils::make_displacement(f, secondPhase); + const mfem::Vector secondField = field_dof_test_utils::make_supported_displacement(f, secondPhase); direction -= secondField; return direction; @@ -93,14 +97,12 @@ namespace prepared_hydrostatic_displacement_test_utils { mfem::Vector residualPlus; mfem::Vector residualMinus; - mean_field::operators::kernels::apply_hydrostatic_equilibrium( - f, *f.domainMapperStateless, rotation, enthalpy, gravityPotential, displacementPlus, bernoulliConstant, - residualPlus + field_dof_test_utils::apply_hydrostatic_reference( + f, rotation, enthalpy, gravityPotential, displacementPlus, bernoulliConstant, residualPlus ); - mean_field::operators::kernels::apply_hydrostatic_equilibrium( - f, *f.domainMapperStateless, rotation, enthalpy, gravityPotential, displacementMinus, bernoulliConstant, - residualMinus + field_dof_test_utils::apply_hydrostatic_reference( + f, rotation, enthalpy, gravityPotential, displacementMinus, bernoulliConstant, residualMinus ); difference = residualPlus; @@ -119,7 +121,7 @@ namespace prepared_hydrostatic_displacement_test_utils { TEST_CASE( "Prepared Hydrostatic Displacement Jacobian Matches Centered Differences", - tags::barotrope &tags::hydro &tags::jacobian &tags::prepared &tags::unit + tags::barotrope_hydrostatic_prepared_jacobian &tags::unit ) { auto args = test_utils::setup_args(); @@ -131,7 +133,7 @@ TEST_CASE( const mfem::Vector gravityPotential = prepared_hydrostatic_displacement_test_utils::make_gravity_potential(f); - const mfem::Vector displacement = gravity_prepared_test_utils::make_displacement(f, 0.73); + const mfem::Vector displacement = field_dof_test_utils::make_supported_displacement(f, 0.73); constexpr double bernoulliConstant = 0.39; @@ -201,7 +203,7 @@ TEST_CASE( TEST_CASE( "Prepared Hydrostatic Displacement Jacobian Reuses And Refreshes Frozen " "Data", - tags::barotrope &tags::hydro &tags::jacobian &tags::prepared &tags::unit + tags::barotrope_hydrostatic_prepared_jacobian &tags::unit ) { auto args = test_utils::setup_args(); @@ -213,7 +215,7 @@ TEST_CASE( const mfem::Vector gravityPotential = prepared_hydrostatic_displacement_test_utils::make_gravity_potential(f, 0.53); - mfem::Vector displacement = gravity_prepared_test_utils::make_displacement(f, 0.42); + mfem::Vector displacement = field_dof_test_utils::make_supported_displacement(f, 0.42); constexpr double bernoulliConstant = 0.36; @@ -284,7 +286,7 @@ TEST_CASE( CHECK(rotationReferenceError < 2.0e-7); CHECK(rotationEffect > 1.0e-8); - displacement = gravity_prepared_test_utils::make_displacement(f, 0.91); + displacement = field_dof_test_utils::make_supported_displacement(f, 0.91); ++dependencies.displacement.revision; diff --git a/tests/operators/prepared_hydrostatic_equilibrium_jacobian.cpp b/tests/operators/prepared_hydrostatic_equilibrium_jacobian.cpp index a8ea0dc..36f3406 100644 --- a/tests/operators/prepared_hydrostatic_equilibrium_jacobian.cpp +++ b/tests/operators/prepared_hydrostatic_equilibrium_jacobian.cpp @@ -34,14 +34,18 @@ namespace prepared_hydrostatic_jacobian_test_utils { const mean_field::fem::FEM &f, const double phase = 0.23 ) { - return gravity_prepared_test_utils::make_deterministic_vector(f.enthalpyFes->GetTrueVSize(), phase); + return field_dof_test_utils::make_deterministic_supported_vector( + *f.enthalpyFes, phase + ); } mfem::Vector make_gravity_potential( const mean_field::fem::FEM &f, const double phase = 0.41 ) { - return gravity_prepared_test_utils::make_deterministic_vector(f.gravityPotentialFes->GetTrueVSize(), phase); + return field_dof_test_utils::make_deterministic_supported_vector( + *f.gravityPotentialFes, phase + ); } mean_field::physics::RigidRotation make_rotation(const double scale = 1.0) { @@ -75,14 +79,12 @@ namespace prepared_hydrostatic_jacobian_test_utils { mfem::Vector residualPlus; mfem::Vector residualMinus; - mean_field::operators::kernels::apply_hydrostatic_equilibrium( - f, *f.domainMapperStateless, rotation, enthalpyPlus, gravityPotentialPlus, displacement, - bernoulliConstantPlus, residualPlus + field_dof_test_utils::apply_hydrostatic_reference( + f, rotation, enthalpyPlus, gravityPotentialPlus, displacement, bernoulliConstantPlus, residualPlus ); - mean_field::operators::kernels::apply_hydrostatic_equilibrium( - f, *f.domainMapperStateless, rotation, enthalpyMinus, gravityPotentialMinus, displacement, - bernoulliConstantMinus, residualMinus + field_dof_test_utils::apply_hydrostatic_reference( + f, rotation, enthalpyMinus, gravityPotentialMinus, displacement, bernoulliConstantMinus, residualMinus ); // The two states are separated by one complete variation: @@ -103,7 +105,7 @@ namespace prepared_hydrostatic_jacobian_test_utils { TEST_CASE( "Prepared Hydrostatic Algebraic Jacobian Matches Centered Residual " "Differences", - tags::barotrope &tags::hydro &tags::jacobian &tags::prepared &tags::unit + tags::barotrope_hydrostatic_prepared_jacobian &tags::unit ) { auto args = test_utils::setup_args(); @@ -115,7 +117,7 @@ TEST_CASE( const mfem::Vector gravityPotential = prepared_hydrostatic_jacobian_test_utils::make_gravity_potential(f); - const mfem::Vector displacement = gravity_prepared_test_utils::make_displacement(f, 0.73); + const mfem::Vector displacement = field_dof_test_utils::make_supported_displacement(f, 0.73); constexpr double bernoulliConstant = 0.39; @@ -251,7 +253,7 @@ TEST_CASE( TEST_CASE( "Prepared Hydrostatic Algebraic Jacobian Reuses And Rebuilds Only With " "Geometry", - tags::barotrope &tags::hydro &tags::jacobian &tags::prepared &tags::unit + tags::barotrope_hydrostatic_prepared_jacobian &tags::unit ) { auto args = test_utils::setup_args(); @@ -263,7 +265,7 @@ TEST_CASE( mfem::Vector gravityPotential = prepared_hydrostatic_jacobian_test_utils::make_gravity_potential(f); - mfem::Vector displacement = gravity_prepared_test_utils::make_displacement(f, 0.42); + mfem::Vector displacement = field_dof_test_utils::make_supported_displacement(f, 0.42); double bernoulliConstant = 0.37; @@ -340,7 +342,7 @@ TEST_CASE( CHECK(preparedOperator.GetAlgebraicJacobianStatistics().preparations == 1); - displacement = gravity_prepared_test_utils::make_displacement(f, 0.73); + displacement = field_dof_test_utils::make_supported_displacement(f, 0.73); ++dependencies.displacement.revision; diff --git a/tests/operators/prepared_mass_normalization.cpp b/tests/operators/prepared_mass_normalization.cpp index 59ba16f..f09b390 100644 --- a/tests/operators/prepared_mass_normalization.cpp +++ b/tests/operators/prepared_mass_normalization.cpp @@ -26,14 +26,25 @@ namespace mass_normalization_test_utils { ); [[nodiscard]] mean_field::operators::MassNormalizationLayout make_layout(const mean_field::fem::FEM &f) { + const auto densityMap = + field_dof_test_utils::make_map(*f.densityFes); + const auto displacementMap = + field_dof_test_utils::make_map(*f.displacementFes); + const auto gravityFluxMap = + field_dof_test_utils::make_map(*f.gravityFluxFes); + const auto gravityPotentialMap = + field_dof_test_utils::make_map(*f.gravityPotentialFes); + const auto enthalpyMap = + field_dof_test_utils::make_map(*f.enthalpyFes); + const std::array valueSizes{ - f.densityFes->GetTrueVSize(), f.displacementFes->GetTrueVSize(), f.gravityFluxFes->GetTrueVSize(), - f.gravityPotentialFes->GetTrueVSize(), f.enthalpyFes->GetTrueVSize(), 1 + densityMap.reduced_size(), displacementMap.reduced_size(), gravityFluxMap.reduced_size(), + gravityPotentialMap.reduced_size(), enthalpyMap.reduced_size(), 1 }; const std::array residualSizes{ - f.gravityFluxFes->GetTrueVSize(), f.gravityPotentialFes->GetTrueVSize(), f.densityFes->GetTrueVSize(), - f.displacementFes->GetTrueVSize(), f.enthalpyFes->GetTrueVSize(), 1 + gravityFluxMap.reduced_size(), gravityPotentialMap.reduced_size(), densityMap.reduced_size(), + displacementMap.reduced_size(), enthalpyMap.reduced_size(), 1 }; return {valueSizes, residualSizes}; @@ -163,10 +174,10 @@ namespace mass_normalization_test_utils { gravityPotential = 0.0; context.Prepare( - {.density = density, - .displacement = displacement, - .gravity_gradient = gravityGradient, - .gravity_potential = gravityPotential}, + {.density = context.GetDensityMap().gather(density), + .displacement = context.GetDisplacementMap().gather(displacement), + .gravity_gradient = context.GetGravityGradientMap().gather(gravityGradient), + .gravity_potential = context.GetGravityPotentialMap().gather(gravityPotential)}, make_gravity_revisions(dependencies, gravityGradientRevision, gravityPotentialRevision) ); } @@ -189,7 +200,7 @@ namespace mass_normalization_test_utils { TEST_CASE( "Prepared Mass Normalization Has The Analytic Affine Volume Scaling", - tags::barotrope &tags::prepared &tags::analytic_comparison + tags::barotrope_mass_normalization_analytic ) { using Operator = mean_field::operators::PreparedMassNormalizationOperator; @@ -265,7 +276,7 @@ TEST_CASE( TEST_CASE( "Prepared Mass Normalization Density Jacobian Matches Centered Difference", - tags::barotrope &tags::prepared &tags::jacobian &tags::accuracy + tags::barotrope_mass_normalization_jacobian &tags::accuracy ) { mean_field::utils::Args args = test_utils::setup_args(); mean_field::fem::FEM f = mean_field::fem::setup_fem(args.mesh_file, args, 0); @@ -285,8 +296,10 @@ TEST_CASE( mean_field::operators::PreparedMassNormalizationOperator massOperator(f, *f.domainMapperStateless, gravityContext); massOperator.Prepare({.targetMass = 1.23}, dependencies); + const mfem::Vector reducedDensityDirection = gravityContext.GetDensityMap().gather(densityDirection); + mfem::Vector analyticAction; - massOperator.ApplyDensityJacobianAction(densityDirection, analyticAction); + massOperator.ApplyDensityJacobianAction(reducedDensityDirection, analyticAction); constexpr double epsilon = 1.0e-3; mfem::Vector densityPlus(density); @@ -312,7 +325,7 @@ TEST_CASE( TEST_CASE( "Prepared Mass Normalization Geometry Jacobian Matches Centered Difference", - tags::barotrope &tags::prepared &tags::jacobian &tags::geometry + tags::barotrope_mass_normalization_jacobian &tags::geometry ) { mean_field::utils::Args args = test_utils::setup_args(); mean_field::fem::FEM f = mean_field::fem::setup_fem(args.mesh_file, args, 0); @@ -332,8 +345,11 @@ TEST_CASE( mean_field::operators::PreparedMassNormalizationOperator massOperator(f, *f.domainMapperStateless, gravityContext); massOperator.Prepare({.targetMass = 1.11}, dependencies); + const mfem::Vector reducedDisplacementDirection = + gravityContext.GetDisplacementMap().gather(displacementDirection); + mfem::Vector analyticAction; - massOperator.ApplyDisplacementJacobianAction(displacementDirection, analyticAction); + massOperator.ApplyDisplacementJacobianAction(reducedDisplacementDirection, analyticAction); constexpr double epsilon = 1.0e-6; mfem::Vector displacementPlus(displacement); @@ -359,7 +375,7 @@ TEST_CASE( TEST_CASE( "Prepared Mass Normalization Selectively Refreshes Its Cached State", - tags::barotrope &tags::prepared &tags::contexts &tags::integration + tags::barotrope_mass_normalization_context &tags::integration ) { mean_field::utils::Args args = test_utils::setup_args(); mean_field::fem::FEM f = mean_field::fem::setup_fem(args.mesh_file, args, 0); @@ -439,7 +455,7 @@ TEST_CASE( TEST_CASE( "Prepared Mass Normalization Complete Action And Coupled Routing Are Exact", - tags::barotrope &tags::prepared &tags::jacobian &tags::mfem_operators + tags::barotrope_mass_normalization_jacobian ) { mean_field::utils::Args args = test_utils::setup_args(); mean_field::fem::FEM f = mean_field::fem::setup_fem(args.mesh_file, args, 0); @@ -460,13 +476,19 @@ TEST_CASE( mean_field::operators::PreparedMassNormalizationOperator massOperator(f, *f.domainMapperStateless, gravityContext); massOperator.Prepare({.targetMass = 1.19}, dependencies); + const mfem::Vector reducedDensityDirection = gravityContext.GetDensityMap().gather(densityDirection); + const mfem::Vector reducedDisplacementDirection = + gravityContext.GetDisplacementMap().gather(displacementDirection); + mfem::Vector densityAction; mfem::Vector displacementAction; mfem::Vector completeAction; - massOperator.ApplyDensityJacobianAction(densityDirection, densityAction); - massOperator.ApplyDisplacementJacobianAction(displacementDirection, displacementAction); - massOperator.ApplyCompleteJacobianAction(densityDirection, displacementDirection, completeAction); + massOperator.ApplyDensityJacobianAction(reducedDensityDirection, densityAction); + massOperator.ApplyDisplacementJacobianAction(reducedDisplacementDirection, displacementAction); + massOperator.ApplyCompleteJacobianAction( + reducedDensityDirection, reducedDisplacementDirection, completeAction + ); CHECK( mass_normalization_test_utils::relative_error(completeAction(0), densityAction(0) + displacementAction(0)) < @@ -476,16 +498,26 @@ TEST_CASE( const auto layout = mass_normalization_test_utils::make_layout(f); mean_field::operators::PreparedMassNormalizationJacobianOperator adapter(layout, massOperator); + CHECK( + layout.size(mass_normalization_test_utils::densityValue) == + gravityContext.GetDensityMap().reduced_size() + ); + CHECK( + layout.size(mass_normalization_test_utils::displacementValue) == + gravityContext.GetDisplacementMap().reduced_size() + ); + CHECK(gravityContext.GetDensityMap().inactive_size() > 0); + mfem::Vector direction(layout.value_offsets().Last()); direction = 0.0; - for (int entry = 0; entry < densityDirection.Size(); ++entry) { - direction(layout.offset(mass_normalization_test_utils::densityValue) + entry) = densityDirection(entry); + for (int entry = 0; entry < reducedDensityDirection.Size(); ++entry) { + direction(layout.offset(mass_normalization_test_utils::densityValue) + entry) = reducedDensityDirection(entry); } - for (int entry = 0; entry < displacementDirection.Size(); ++entry) { + for (int entry = 0; entry < reducedDisplacementDirection.Size(); ++entry) { direction(layout.offset(mass_normalization_test_utils::displacementValue) + entry) = - displacementDirection(entry); + reducedDisplacementDirection(entry); } mfem::Vector coupledAction; @@ -505,4 +537,4 @@ TEST_CASE( CHECK(&massOperator.GetFEM() == &f); CHECK(&massOperator.GetGravityContext() == &gravityContext); CHECK(adapter.GetLayout().residual_offsets().Last() == layout.residual_offsets().Last()); -} \ No newline at end of file +} diff --git a/tests/operators/prepared_rotation_displacement_force.cpp b/tests/operators/prepared_rotation_displacement_force.cpp index 4169844..5c15a34 100644 --- a/tests/operators/prepared_rotation_displacement_force.cpp +++ b/tests/operators/prepared_rotation_displacement_force.cpp @@ -58,14 +58,25 @@ namespace rotational_displacement_force_test_utils { ); [[nodiscard]] mean_field::operators::RotationalDisplacementForceLayout make_layout(const mean_field::fem::FEM &f) { + using DomainSchema = gravity_prepared_test_utils::DomainSchema; + + const auto densityMap = gravity_prepared_test_utils::make_field_map(f); + const auto displacementMap = gravity_prepared_test_utils::make_field_map(f); + const auto gravityFluxMap = + mean_field::field::make_field_dof_map(*f.gravityFluxFes); + const auto gravityPotentialMap = + mean_field::field::make_field_dof_map(*f.gravityPotentialFes); + const auto enthalpyMap = + mean_field::field::make_field_dof_map(*f.enthalpyFes); + const std::array valueSizes{ - f.densityFes->GetTrueVSize(), f.displacementFes->GetTrueVSize(), f.gravityFluxFes->GetTrueVSize(), - f.gravityPotentialFes->GetTrueVSize(), f.enthalpyFes->GetTrueVSize(), 1 + densityMap.reduced_size(), displacementMap.reduced_size(), gravityFluxMap.reduced_size(), + gravityPotentialMap.reduced_size(), enthalpyMap.reduced_size(), 1 }; const std::array residualSizes{ - f.gravityFluxFes->GetTrueVSize(), f.gravityPotentialFes->GetTrueVSize(), f.densityFes->GetTrueVSize(), - f.displacementFes->GetTrueVSize(), f.enthalpyFes->GetTrueVSize(), 1 + gravityFluxMap.reduced_size(), gravityPotentialMap.reduced_size(), densityMap.reduced_size(), + displacementMap.reduced_size(), enthalpyMap.reduced_size(), 1 }; return {valueSizes, residualSizes}; @@ -275,7 +286,7 @@ namespace rotational_displacement_force_test_utils { TEST_CASE( "Rotational Displacement Force Query Includes Density Test And Linear " "Position", - tags::centrifugal &tags::quadrature &tags::unit + tags::rotation_prepared_unit ) { using DisplacementField = mean_field::field::Field; @@ -300,7 +311,7 @@ TEST_CASE( TEST_CASE( "Rotational Displacement Force Uses Negative Rotation-Potential " "Gradient And Excludes Vacuum", - tags::centrifugal &tags::kernels &tags::integration &tags::accuracy + tags::rotation_kernel_accuracy ) { mean_field::utils::Args args = test_utils::setup_args(); @@ -368,7 +379,7 @@ TEST_CASE( TEST_CASE( "Prepared Rotational Displacement Force Reprepares Selectively", - tags::centrifugal &tags::prepared &tags::integration + tags::rotation_prepared ) { mean_field::utils::Args args = test_utils::setup_args(); @@ -376,9 +387,9 @@ TEST_CASE( REQUIRE(f.okay()); - mfem::Vector density = rotational_displacement_force_test_utils::make_density(f, 0.37); + mfem::Vector densityTrue = rotational_displacement_force_test_utils::make_density(f, 0.37); - const mfem::Vector displacement = gravity_prepared_test_utils::make_displacement(f, 0.53); + const mfem::Vector displacementTrue = gravity_prepared_test_utils::make_displacement(f, 0.53); mean_field::physics::RigidRotation rotation = rotational_displacement_force_test_utils::make_rotation(0.81); @@ -386,6 +397,10 @@ TEST_CASE( mean_field::operators::PreparedRotationalDisplacementForceOperator preparedOperator(f, *f.domainMapperStateless); + const auto &context = preparedOperator.GetContext(); + mfem::Vector density = context.GetDensityMap().gather(densityTrue); + const mfem::Vector displacement = context.GetDisplacementMap().gather(displacementTrue); + const auto initialReport = preparedOperator.Prepare({.density = density, .displacement = displacement}, dependencies, rotation); @@ -400,19 +415,22 @@ TEST_CASE( preparedOperator.BuildResidual(preparedResidual); mean_field::operators::kernels::apply_rotational_displacement_force_residual( - f, *f.domainMapperStateless, rotation, density, displacement, kernelResidual + f, *f.domainMapperStateless, rotation, densityTrue, displacementTrue, kernelResidual ); + const mfem::Vector kernelResidualReduced = context.GetDisplacementMap().gather(kernelResidual); + CHECK( rotational_displacement_force_test_utils::relative_difference( - preparedResidual, kernelResidual, f.mesh->GetComm() + preparedResidual, kernelResidualReduced, f.mesh->GetComm() ) < 2.0e-12 ); CHECK_FALSE(preparedOperator.Prepare({.density = density, .displacement = displacement}, dependencies, rotation) .DidAnyWork()); - density = rotational_displacement_force_test_utils::make_density(f, 0.79); + densityTrue = rotational_displacement_force_test_utils::make_density(f, 0.79); + density = context.GetDensityMap().gather(densityTrue); ++dependencies.density.revision; @@ -438,7 +456,7 @@ TEST_CASE( TEST_CASE( "Rotational Displacement Force Jacobian Matches Both Columns And " "Centered Differences", - tags::centrifugal &tags::prepared &tags::jacobian &tags::accuracy + tags::rotation_prepared_jacobian_accuracy ) { mean_field::utils::Args args = test_utils::setup_args(); @@ -446,18 +464,25 @@ TEST_CASE( REQUIRE(f.okay()); - const mfem::Vector density = rotational_displacement_force_test_utils::make_density(f, 0.43); + const mfem::Vector densityTrue = rotational_displacement_force_test_utils::make_density(f, 0.43); - const mfem::Vector densityDirection = rotational_displacement_force_test_utils::make_density_direction(f, 0.59); + const mfem::Vector densityDirectionTrue = rotational_displacement_force_test_utils::make_density_direction(f, 0.59); - const mfem::Vector displacement = gravity_prepared_test_utils::make_displacement(f, 0.61); + const mfem::Vector displacementTrue = gravity_prepared_test_utils::make_displacement(f, 0.61); - const mfem::Vector displacementDirection = rotational_displacement_force_test_utils::make_displacement_direction(f); + const mfem::Vector displacementDirectionTrue = + rotational_displacement_force_test_utils::make_displacement_direction(f); const mean_field::physics::RigidRotation rotation = rotational_displacement_force_test_utils::make_rotation(0.93); mean_field::operators::PreparedRotationalDisplacementForceOperator preparedOperator(f, *f.domainMapperStateless); + const auto &context = preparedOperator.GetContext(); + const mfem::Vector density = context.GetDensityMap().gather(densityTrue); + const mfem::Vector densityDirection = context.GetDensityMap().gather(densityDirectionTrue); + const mfem::Vector displacement = context.GetDisplacementMap().gather(displacementTrue); + const mfem::Vector displacementDirection = context.GetDisplacementMap().gather(displacementDirectionTrue); + preparedOperator.Prepare( {.density = density, .displacement = displacement}, rotational_displacement_force_test_utils::make_dependencies(), rotation @@ -482,26 +507,30 @@ TEST_CASE( ) < 2.0e-12 ); - mfem::Vector zeroDensity(densityDirection.Size()); - mfem::Vector zeroDisplacement(displacementDirection.Size()); - zeroDensity = 0.0; - zeroDisplacement = 0.0; + mfem::Vector zeroDensityTrue(densityDirectionTrue.Size()); + mfem::Vector zeroDisplacementTrue(displacementDirectionTrue.Size()); + zeroDensityTrue = 0.0; + zeroDisplacementTrue = 0.0; - constexpr double step = 1.0e-5; + constexpr double step = 1.0e-5; - const mfem::Vector densityDifference = rotational_displacement_force_test_utils::centered_difference( - f, rotation, density, densityDirection, displacement, zeroDisplacement, step + const mfem::Vector densityDifferenceTrue = rotational_displacement_force_test_utils::centered_difference( + f, rotation, densityTrue, densityDirectionTrue, displacementTrue, zeroDisplacementTrue, step ); - const mfem::Vector displacementDifference = rotational_displacement_force_test_utils::centered_difference( - f, rotation, density, zeroDensity, displacement, displacementDirection, step + const mfem::Vector displacementDifferenceTrue = rotational_displacement_force_test_utils::centered_difference( + f, rotation, densityTrue, zeroDensityTrue, displacementTrue, displacementDirectionTrue, step ); - const mfem::Vector completeDifference = rotational_displacement_force_test_utils::centered_difference( - f, rotation, density, densityDirection, displacement, displacementDirection, step + const mfem::Vector completeDifferenceTrue = rotational_displacement_force_test_utils::centered_difference( + f, rotation, densityTrue, densityDirectionTrue, displacementTrue, displacementDirectionTrue, step ); - const double densityError = rotational_displacement_force_test_utils::relative_difference( + const mfem::Vector densityDifference = context.GetDisplacementMap().gather(densityDifferenceTrue); + const mfem::Vector displacementDifference = context.GetDisplacementMap().gather(displacementDifferenceTrue); + const mfem::Vector completeDifference = context.GetDisplacementMap().gather(completeDifferenceTrue); + + const double densityError = rotational_displacement_force_test_utils::relative_difference( densityAction, densityDifference, f.mesh->GetComm() ); @@ -524,7 +553,7 @@ TEST_CASE( TEST_CASE( "Prepared Rotational Displacement Force MFEM Adapter Routes Only R-d", - tags::centrifugal &tags::prepared &tags::mfem_operators &tags::unit + tags::rotation_prepared_unit ) { mean_field::utils::Args args = test_utils::setup_args(); @@ -532,18 +561,25 @@ TEST_CASE( REQUIRE(f.okay()); - const mfem::Vector density = rotational_displacement_force_test_utils::make_density(f, 0.47); + const mfem::Vector densityTrue = rotational_displacement_force_test_utils::make_density(f, 0.47); - const mfem::Vector densityDirection = rotational_displacement_force_test_utils::make_density_direction(f, 0.63); + const mfem::Vector densityDirectionTrue = rotational_displacement_force_test_utils::make_density_direction(f, 0.63); - const mfem::Vector displacement = gravity_prepared_test_utils::make_displacement(f, 0.57); + const mfem::Vector displacementTrue = gravity_prepared_test_utils::make_displacement(f, 0.57); - const mfem::Vector displacementDirection = rotational_displacement_force_test_utils::make_displacement_direction(f); + const mfem::Vector displacementDirectionTrue = + rotational_displacement_force_test_utils::make_displacement_direction(f); const mean_field::physics::RigidRotation rotation = rotational_displacement_force_test_utils::make_rotation(0.87); mean_field::operators::PreparedRotationalDisplacementForceOperator preparedOperator(f, *f.domainMapperStateless); + const auto &context = preparedOperator.GetContext(); + const mfem::Vector density = context.GetDensityMap().gather(densityTrue); + const mfem::Vector densityDirection = context.GetDensityMap().gather(densityDirectionTrue); + const mfem::Vector displacement = context.GetDisplacementMap().gather(displacementTrue); + const mfem::Vector displacementDirection = context.GetDisplacementMap().gather(displacementDirectionTrue); + preparedOperator.Prepare( {.density = density, .displacement = displacement}, rotational_displacement_force_test_utils::make_dependencies(), rotation diff --git a/tests/operators/prepared_rotation_displacement_force_affine_deformation.cpp b/tests/operators/prepared_rotation_displacement_force_affine_deformation.cpp index 07ffb00..50186a9 100644 --- a/tests/operators/prepared_rotation_displacement_force_affine_deformation.cpp +++ b/tests/operators/prepared_rotation_displacement_force_affine_deformation.cpp @@ -152,7 +152,7 @@ namespace rotational_displacement_force_affine_deformation_test_utils { TEST_CASE( "Rotational Displacement Force Matches A Nontrivially Deformed " "Homogeneous Ellipsoid", - tags::centrifugal &tags::analytic_comparison &tags::accuracy &tags::geometry + tags::rotation_analytic_accuracy_geometry ) { mean_field::utils::Args args = test_utils::setup_args(); diff --git a/tests/operators/prepared_rotation_displacement_force_analytic.cpp b/tests/operators/prepared_rotation_displacement_force_analytic.cpp index afa238c..5fb1318 100644 --- a/tests/operators/prepared_rotation_displacement_force_analytic.cpp +++ b/tests/operators/prepared_rotation_displacement_force_analytic.cpp @@ -67,7 +67,7 @@ namespace rotational_displacement_force_analytic_test_utils { TEST_CASE( "Rigid Rotation Gradient And Hessian Action Match The Analytic " "Potential", - tags::centrifugal &tags::unit &tags::accuracy + tags::rotation_analytic_unit ) { mfem::Vector angularVelocity(3); angularVelocity(0) = 0.23; @@ -130,7 +130,7 @@ TEST_CASE( TEST_CASE( "Rotational Displacement Force Reproduces The Homogeneous Sphere " "Rotational Virial", - tags::centrifugal &tags::analytic_comparison &tags::accuracy + tags::rotation_analytic_accuracy ) { mean_field::utils::Args args = test_utils::setup_args(); @@ -187,7 +187,7 @@ TEST_CASE( TEST_CASE( "Rotational Displacement Force Matches The Analytic Off-Axis " "Resultant", - tags::centrifugal &tags::analytic_comparison &tags::accuracy + tags::rotation_analytic_accuracy ) { mean_field::utils::Args args = test_utils::setup_args(); diff --git a/tests/operators/prepared_stellar_equilibrium.cpp b/tests/operators/prepared_stellar_equilibrium.cpp index 3559b07..1403c26 100644 --- a/tests/operators/prepared_stellar_equilibrium.cpp +++ b/tests/operators/prepared_stellar_equilibrium.cpp @@ -153,23 +153,16 @@ namespace stellar_equilibrium_test_utils { return map.gather(fullEnthalpy); } - [[nodiscard]] mfem::Vector pack_full_gravity_state( - const mean_field::fem::FEM &f, - const mfem::Vector &fullDensity, + [[nodiscard]] mfem::Vector pack_gravity_state( + const mfem::Vector &density, const mfem::Vector &displacement, const mfem::Vector &gravityGradient, const mfem::Vector &gravityPotential ) { const std::array blockSizes{ - f.densityFes->GetTrueVSize(), f.displacementFes->GetTrueVSize(), f.gravityFluxFes->GetTrueVSize(), - f.gravityPotentialFes->GetTrueVSize() + density.Size(), displacement.Size(), gravityGradient.Size(), gravityPotential.Size() }; - MFEM_VERIFY(fullDensity.Size() == blockSizes[0], "Full density has the wrong gravity-state size."); - MFEM_VERIFY(displacement.Size() == blockSizes[1], "Displacement has the wrong gravity-state size."); - MFEM_VERIFY(gravityGradient.Size() == blockSizes[2], "Gravity gradient has the wrong gravity-state size."); - MFEM_VERIFY(gravityPotential.Size() == blockSizes[3], "Gravity potential has the wrong gravity-state size."); - const std::array offsets{ 0, blockSizes[0], blockSizes[0] + blockSizes[1], blockSizes[0] + blockSizes[1] + blockSizes[2], blockSizes[0] + blockSizes[1] + blockSizes[2] + blockSizes[3] @@ -178,7 +171,7 @@ namespace stellar_equilibrium_test_utils { mfem::Vector packed(offsets[4]); const std::array blocks{ - &fullDensity, &displacement, &gravityGradient, &gravityPotential + &density, &displacement, &gravityGradient, &gravityPotential }; for (int block = 0; block < 4; ++block) { @@ -494,27 +487,24 @@ namespace stellar_equilibrium_test_utils { const mfem::Vector &state ) { const mean_field::operators::StellarEquilibriumLayout &layout = stellarOperator.GetLayout(); - const FieldMaps maps(f); - const mfem::Vector reducedDensity = const_value_view(state, layout, densityValue); const mfem::Vector displacement = const_value_view(state, layout, displacementValue); const mfem::Vector gravityGradient = const_value_view(state, layout, gravityGradientValue); const mfem::Vector gravityPotential = const_value_view(state, layout, gravityPotentialValue); - const mfem::Vector fullDensity = maps.density.scatter(reducedDensity); const mfem::Vector gravityState = - pack_full_gravity_state(f, fullDensity, displacement, gravityGradient, gravityPotential); + pack_gravity_state(reducedDensity, displacement, gravityGradient, gravityPotential); mfem::Vector gravity; mfem::Vector closure; mfem::Vector displacementResidualValue; - mfem::Vector fullHydrostatic; + mfem::Vector hydrostatic; mfem::Vector mass; stellarOperator.GetGravityOperator().Mult(gravityState, gravity); stellarOperator.GetBarotropicClosureOperator().BuildResidual(closure); stellarOperator.GetDisplacementOperator().BuildResidual(displacementResidualValue); - stellarOperator.GetHydrostaticOperator().BuildResidual(fullHydrostatic); + stellarOperator.GetHydrostaticOperator().BuildResidual(hydrostatic); stellarOperator.GetMassNormalizationOperator().BuildResidual(mass); mfem::Vector result(layout.residual_offsets().Last()); @@ -537,10 +527,7 @@ namespace stellar_equilibrium_test_utils { residual_view(result, layout, displacementResidual) = displacementResidualValue; - { - mfem::Vector reducedHydrostatic = residual_view(result, layout, enthalpyResidual); - maps.enthalpy.gather(fullHydrostatic, reducedHydrostatic); - } + residual_view(result, layout, enthalpyResidual) = hydrostatic; residual_view(result, layout, massResidual) = mass; @@ -549,11 +536,9 @@ namespace stellar_equilibrium_test_utils { [[nodiscard]] mfem::Vector explicit_jacobian_action( const mean_field::operators::PreparedStellarEquilibriumOperator &stellarOperator, - const mean_field::fem::FEM &f, const mfem::Vector &direction ) { const mean_field::operators::StellarEquilibriumLayout &layout = stellarOperator.GetLayout(); - const FieldMaps maps(f); const mfem::Vector reducedDensityDirection = const_value_view(direction, layout, densityValue); const mfem::Vector displacementDirection = const_value_view(direction, layout, displacementValue); @@ -562,17 +547,14 @@ namespace stellar_equilibrium_test_utils { const mfem::Vector reducedEnthalpyDirection = const_value_view(direction, layout, enthalpyValue); const mfem::Vector bernoulliDirection = const_value_view(direction, layout, bernoulliValue); - const mfem::Vector fullDensityDirection = maps.density.scatter(reducedDensityDirection); - const mfem::Vector fullEnthalpyDirection = maps.enthalpy.scatter(reducedEnthalpyDirection); - - const mfem::Vector gravityDirection = pack_full_gravity_state( - f, fullDensityDirection, displacementDirection, gravityGradientDirection, gravityPotentialDirection + const mfem::Vector gravityDirection = pack_gravity_state( + reducedDensityDirection, displacementDirection, gravityGradientDirection, gravityPotentialDirection ); mfem::Vector gravityAction; mfem::Vector closureAction; mfem::Vector displacementAction; - mfem::Vector fullHydrostaticAction; + mfem::Vector hydrostaticAction; mfem::Vector massAction; stellarOperator.GetGravityJacobianOperator().Mult(gravityDirection, gravityAction); @@ -582,17 +564,17 @@ namespace stellar_equilibrium_test_utils { ); stellarOperator.GetDisplacementOperator().ApplyCompleteJacobianAction( - fullDensityDirection, displacementDirection, gravityGradientDirection, fullEnthalpyDirection, + reducedDensityDirection, displacementDirection, gravityGradientDirection, reducedEnthalpyDirection, displacementAction ); stellarOperator.GetHydrostaticOperator().ApplyCompleteJacobianAction( - fullEnthalpyDirection, gravityPotentialDirection, bernoulliDirection(0), displacementDirection, - fullHydrostaticAction + reducedEnthalpyDirection, gravityPotentialDirection, bernoulliDirection(0), displacementDirection, + hydrostaticAction ); stellarOperator.GetMassNormalizationOperator().ApplyCompleteJacobianAction( - fullDensityDirection, displacementDirection, massAction + reducedDensityDirection, displacementDirection, massAction ); mfem::Vector result(layout.residual_offsets().Last()); @@ -615,10 +597,7 @@ namespace stellar_equilibrium_test_utils { residual_view(result, layout, displacementResidual) = displacementAction; - { - mfem::Vector reducedHydrostaticAction = residual_view(result, layout, enthalpyResidual); - maps.enthalpy.gather(fullHydrostaticAction, reducedHydrostaticAction); - } + residual_view(result, layout, enthalpyResidual) = hydrostaticAction; residual_view(result, layout, massResidual) = massAction; @@ -719,7 +698,7 @@ TEST_CASE( TEST_CASE( "Prepared Stellar Equilibrium Jacobian Is The Exact Restricted Full Child Jacobian", - tags::barotrope &tags::prepared &tags::field &tags::jacobian &tags::integration &tags::accuracy + tags::barotrope_prepared_jacobian_accuracy &tags::field ) { mean_field::utils::Args args = test_utils::setup_args(); mean_field::fem::FEM f = mean_field::fem::setup_fem(args.mesh_file, args, 0); @@ -742,7 +721,7 @@ TEST_CASE( stellarOperator.Mult(direction, rootAction); const mfem::Vector explicitAction = - stellar_equilibrium_test_utils::explicit_jacobian_action(stellarOperator, f, direction); + stellar_equilibrium_test_utils::explicit_jacobian_action(stellarOperator, direction); const double difference = stellar_equilibrium_test_utils::relative_difference(rootAction, explicitAction, f.mesh->GetComm()); diff --git a/tests/physics/gravity.cpp b/tests/physics/gravity.cpp index 083b46a..a255001 100644 --- a/tests/physics/gravity.cpp +++ b/tests/physics/gravity.cpp @@ -1077,7 +1077,7 @@ namespace { TEST_CASE( "Uniform Potential Matches Analytic", - tags::gravity &tags::analytic_comparison &tags::initialization + tags::gravity_analytic_initialization ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -1161,7 +1161,7 @@ TEST_CASE( TEST_CASE( "Parabolic Density Virial Self-Consistency", - tags::gravity &tags::analytic_comparison &tags::initialization + tags::gravity_analytic_initialization ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -1218,7 +1218,7 @@ TEST_CASE( TEST_CASE( "Rational Density Virial Self-Consistency", - tags::gravity &tags::self_consistency &tags::initialization + tags::gravity_consistency_initialization ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -1278,7 +1278,7 @@ TEST_CASE( TEST_CASE( "Homogeneous Ellipsoid Analytic Gravity", - tags::gravity &tags::analytic_comparison &tags::initialization + tags::gravity_analytic_initialization ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -1482,7 +1482,7 @@ TEST_CASE( TEST_CASE( "Deformed Rational Density Virial Self-Consistency", - tags::gravity &tags::self_consistency &tags::initialization + tags::gravity_consistency_initialization ) { auto args = test_utils::setup_args(); fem::FEM f = fem::setup_fem(args.mesh_file, args, 0); @@ -1565,7 +1565,7 @@ TEST_CASE( TEST_CASE( "New Gravity Potential Matches Uniform Sphere Analytic", - tags::gravity &tags::analytic_comparison &tags::initialization &tags::integration + tags::gravity_analytic ) { auto args = test_utils::setup_args(); args.p.rtol = 1.0e-13; @@ -1675,7 +1675,7 @@ TEST_CASE( TEST_CASE( "New Gravity Potential Matches Legacy Solver On Homogeneous Ellipsoid", - tags::gravity &tags::analytic_comparison &tags::initialization &tags::integration &tags::legacy_comparison + tags::gravity_legacy ) { auto args = test_utils::setup_args(); args.p.rtol = 1.0e-13; @@ -1739,13 +1739,22 @@ TEST_CASE( constexpr auto gravity_poisson_residual_block = utils::blocks::get_residual_block(utils::blocks::gravity_field.poisson_term); + using DomainSchema = utils::domain::CoreEnvelopeVacuumDomainSchema; + const field::FieldDofMap density_map = field::make_field_dof_map(*f.densityFes); + const field::FieldDofMap displacement_map = + field::make_field_dof_map(*f.displacementFes); + const field::FieldDofMap gravity_flux_map = + field::make_field_dof_map(*f.gravityFluxFes); + const field::FieldDofMap gravity_potential_map = + field::make_field_dof_map(*f.gravityPotentialFes); + const std::array value_sizes{ - f.densityFes->GetTrueVSize(), f.displacementFes->GetTrueVSize(), f.gravityFluxFes->GetTrueVSize(), - f.gravityPotentialFes->GetTrueVSize() + density_map.reduced_size(), displacement_map.reduced_size(), gravity_flux_map.reduced_size(), + gravity_potential_map.reduced_size() }; const std::array residual_sizes{ - f.gravityFluxFes->GetTrueVSize(), f.gravityPotentialFes->GetTrueVSize() + gravity_flux_map.reduced_size(), gravity_potential_map.reduced_size() }; const utils::blocks::form_layout gravity_layout(value_sizes, residual_sizes); @@ -1760,6 +1769,8 @@ TEST_CASE( REQUIRE(density_true.Size() == f.densityFes->GetTrueVSize()); REQUIRE(displacement_true.Size() == f.displacementFes->GetTrueVSize()); + const mfem::Vector reduced_density = density_map.gather(density_true); + const mfem::Vector reduced_displacement = displacement_map.gather(displacement_true); operators::context::gravity_field::GravityFieldLinearizationContext residual_linearization_context( f, *f.domainMapperStateless @@ -1779,24 +1790,24 @@ TEST_CASE( ); operators::ReducedGravityFieldOperator new_reduced_operator( - residual_gravity_operator, residual_geometry_context, displacement_true + residual_gravity_operator, residual_geometry_context, reduced_displacement ); REQUIRE(new_reduced_operator.Width() == gravity_system_size); REQUIRE(new_reduced_operator.Height() == gravity_system_size); mfem::Vector new_right_hand_side; - new_reduced_operator.BuildRightHandSide(density_true, new_right_hand_side); + new_reduced_operator.BuildRightHandSide(reduced_density, new_right_hand_side); REQUIRE(new_right_hand_side.Size() == gravity_system_size); - REQUIRE(f.gravityContext.source_form->Width() == density_true.Size()); + REQUIRE(f.gravityContext.source_form->Width() == reduced_density.Size()); REQUIRE(f.gravityContext.source_form->Height() == gravity_poisson_size); mfem::Vector legacy_source_action(f.gravityContext.source_form->Height()); legacy_source_action = 0.0; - f.gravityContext.source_form->Mult(density_true, legacy_source_action); + f.gravityContext.source_form->Mult(reduced_density, legacy_source_action); REQUIRE(legacy_source_action.Size() == gravity_poisson_size); @@ -2057,7 +2068,7 @@ TEST_CASE( TEST_CASE( "New Gravity Potential Deformed Rational Density Virial Self-Consistency", - tags::gravity &tags::self_consistency &tags::initialization &tags::integration + tags::gravity_consistency ) { auto args = test_utils::setup_args(); args.p.rtol = 1.0e-13; @@ -2183,7 +2194,7 @@ TEST_CASE( TEST_CASE( "New Gravity Potential Resolves Exterior Monopole By Compactification " "Shell", - tags::gravity &tags::analytic_comparison &tags::initialization &tags::integration + tags::gravity_analytic ) { auto args = test_utils::setup_args(); args.p.rtol = 1.0e-13; @@ -2305,7 +2316,7 @@ TEST_CASE( TEST_CASE( "New Exterior Monopole Error Is Separated From Finite Element Projection " "Floor", - tags::gravity &tags::analytic_comparison &tags::initialization &tags::integration &tags::accuracy + tags::gravity_analytic_accuracy ) { auto args = test_utils::setup_args(); args.p.rtol = 1.0e-13; @@ -2466,7 +2477,7 @@ TEST_CASE( TEST_CASE( "New Gravity Potential Matches Analytic Interior Potential For A Deformed " "Homogeneous Star", - tags::gravity &tags::analytic_comparison &tags::initialization &tags::integration &tags::accuracy + tags::gravity_analytic_accuracy ) { auto args = test_utils::setup_args(); args.p.rtol = 1.0e-13; @@ -2807,7 +2818,7 @@ static void evaluate_ferrers_n1_gradient( TEST_CASE( "New Gravity Potential Matches Analytic Ferrers Ellipsoid", - tags::gravity &tags::analytic_comparison &tags::initialization &tags::integration &tags::accuracy + tags::gravity_analytic_accuracy ) { auto args = test_utils::setup_args(); diff --git a/tests/physics/gravity_monopole_accuracy.cpp b/tests/physics/gravity_monopole_accuracy.cpp index 68635de..7679733 100644 --- a/tests/physics/gravity_monopole_accuracy.cpp +++ b/tests/physics/gravity_monopole_accuracy.cpp @@ -458,9 +458,10 @@ namespace { const mfem::Vector projection_rhs = assemble_monopole_projection_rhs(f, displacement, mass, stellar_radius); mean_field::operators::PreparedMappedHDivMassOperator mass_operator(f, *f.domainMapperStateless); - mass_operator.Prepare(displacement_true); + mass_operator.Prepare(mass_operator.GetDisplacementMap().gather(displacement_true)); + const mfem::Vector reduced_projection_rhs = mass_operator.GetFluxMap().gather(projection_rhs); - mfem::Vector projected_gradient(f.gravityFluxFes->GetTrueVSize()); + mfem::Vector projected_gradient(mass_operator.Width()); projected_gradient = 0.0; mfem::CGSolver solver(f.gravityFluxFes->GetComm()); @@ -469,13 +470,13 @@ namespace { solver.SetAbsTol(1.0e-12); solver.SetMaxIter(2000); solver.SetPrintLevel(0); - solver.Mult(projection_rhs, projected_gradient); + solver.Mult(reduced_projection_rhs, projected_gradient); mfem::Vector projection_residual; mass_operator.Mult(projected_gradient, projection_residual); - projection_residual -= projection_rhs; + projection_residual -= reduced_projection_rhs; - const double source_norm = global_norm(projection_rhs, f.gravityFluxFes->GetComm()); + const double source_norm = global_norm(reduced_projection_rhs, f.gravityFluxFes->GetComm()); const double residual_norm = global_norm(projection_residual, f.gravityFluxFes->GetComm()); const double relative_residual = residual_norm / std::max(source_norm, std::numeric_limits::epsilon()); @@ -487,7 +488,7 @@ namespace { REQUIRE(std::isfinite(relative_residual)); REQUIRE(relative_residual < 1.0e-8); - return projected_gradient; + return mass_operator.GetFluxMap().scatter(projected_gradient); } double mapped_hdiv_relative_error( @@ -502,7 +503,7 @@ namespace { displacement.GetTrueDofs(displacement_true); mean_field::operators::PreparedMappedHDivMassOperator mass_operator(f, *f.domainMapperStateless); - mass_operator.Prepare(displacement_true); + mass_operator.Prepare(mass_operator.GetDisplacementMap().gather(displacement_true)); mfem::Vector difference(computed); difference -= reference; @@ -510,13 +511,15 @@ namespace { mfem::Vector difference_action; mfem::Vector reference_action; - mass_operator.Mult(difference, difference_action); - mass_operator.Mult(reference, reference_action); + const mfem::Vector reduced_difference = mass_operator.GetFluxMap().gather(difference); + const mfem::Vector reduced_reference = mass_operator.GetFluxMap().gather(reference); + mass_operator.Mult(reduced_difference, difference_action); + mass_operator.Mult(reduced_reference, reference_action); MPI_Comm communicator = f.gravityFluxFes->GetComm(); - const double difference_energy = global_dot(difference, difference_action, communicator); - const double reference_energy = global_dot(reference, reference_action, communicator); + const double difference_energy = global_dot(reduced_difference, difference_action, communicator); + const double reference_energy = global_dot(reduced_reference, reference_action, communicator); REQUIRE(difference_energy >= -1.0e-12 * std::abs(reference_energy)); REQUIRE(reference_energy > 0.0); @@ -853,7 +856,7 @@ namespace { TEST_CASE( "New Gravity Monopole Accuracy And Projection Floor", - tags::gravity &tags::analytic_comparison &tags::initialization &tags::integration &tags::accuracy + tags::gravity_analytic_accuracy ) { auto args = test_utils::setup_args(); args.p.rtol = 1.0e-13; @@ -1022,7 +1025,7 @@ TEST_CASE( TEST_CASE( "New Gravity Virial Consistency Across Volume Preserving Deformation", - tags::gravity &tags::self_consistency &tags::initialization &tags::integration &tags::accuracy + tags::gravity_consistency &tags::accuracy ) { auto args = test_utils::setup_args(); args.p.rtol = 1.0e-13; diff --git a/tests/test_helpers.cppm b/tests/test_helpers.cppm index 1c5ede1..9609d05 100644 --- a/tests/test_helpers.cppm +++ b/tests/test_helpers.cppm @@ -2,6 +2,7 @@ module; #include #include #include +#include #include #include @@ -86,6 +87,27 @@ export namespace test_utils { } // namespace test_utils export namespace gravity_prepared_test_utils { + using DomainSchema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema; + + template inline mean_field::field::FieldDofMap make_field_map(const mean_field::fem::FEM &f) { + if constexpr (std::same_as) { + return mean_field::field::make_field_dof_map(*f.densityFes); + } else if constexpr (std::same_as) { + return mean_field::field::make_field_dof_map(*f.displacementFes); + } else { + static_assert(std::same_as); + return mean_field::field::make_field_dof_map(*f.gravityFluxFes); + } + } + + template + inline mfem::Vector gather_field( + const mean_field::fem::FEM &f, + const mfem::Vector &true_vector + ) { + return make_field_map(f).gather(true_vector); + } + inline mfem::Vector make_deterministic_vector( const int size, const double phase = 0.0 @@ -207,59 +229,193 @@ export namespace gravity_prepared_test_utils { } } // namespace gravity_prepared_test_utils +export namespace field_dof_test_utils { + using DomainSchema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema; + + template + inline mean_field::field::FieldDofMap make_map(const mfem::ParFiniteElementSpace &finiteElementSpace) { + return mean_field::field::make_field_dof_map(finiteElementSpace); + } + + template + inline mfem::Vector make_deterministic_supported_vector( + const mfem::ParFiniteElementSpace &finiteElementSpace, + const double phase + ) { + const mean_field::field::FieldDofMap map = make_map(finiteElementSpace); + const mfem::Vector full = + gravity_prepared_test_utils::make_deterministic_vector(map.full_size(), phase); + return map.gather(full); + } + + inline mfem::Vector make_supported_displacement( + const mean_field::fem::FEM &f, + const double phase + ) { + const mean_field::field::FieldDofMap map = + make_map(*f.displacementFes); + return map.gather(gravity_prepared_test_utils::make_displacement(f, phase)); + } + + inline void apply_hydrostatic_reference( + const mean_field::fem::FEM &f, + const mean_field::physics::RigidRotation &rotation, + const mfem::Vector &enthalpy, + const mfem::Vector &gravityPotential, + const mfem::Vector &displacement, + const double bernoulliConstant, + mfem::Vector &residual + ) { + const mean_field::field::FieldDofMap enthalpyMap = + make_map(*f.enthalpyFes); + const mean_field::field::FieldDofMap gravityPotentialMap = + make_map(*f.gravityPotentialFes); + const mean_field::field::FieldDofMap displacementMap = + make_map(*f.displacementFes); + + mfem::Vector enthalpyTrue(enthalpyMap.full_size()); + mfem::Vector gravityPotentialTrue(gravityPotentialMap.full_size()); + mfem::Vector displacementTrue(displacementMap.full_size()); + mfem::Vector residualTrue; + + enthalpyMap.scatter(enthalpy, enthalpyTrue); + gravityPotentialMap.scatter(gravityPotential, gravityPotentialTrue); + displacementMap.scatter(displacement, displacementTrue); + + mean_field::operators::kernels::apply_hydrostatic_equilibrium( + f, *f.domainMapperStateless, rotation, enthalpyTrue, gravityPotentialTrue, displacementTrue, + bernoulliConstant, residualTrue + ); + + residual.SetSize(enthalpyMap.reduced_size()); + enthalpyMap.gather(residualTrue, residual); + } +} // namespace field_dof_test_utils + export namespace tags { - inline constexpr auto geometry = make_tag("geometry"); - inline constexpr auto physics = make_tag("physics"); - inline constexpr auto unit = make_tag("unit"); - inline constexpr auto mesh = make_tag("mesh"); - inline constexpr auto integration = make_tag("integration"); - inline constexpr auto solver = make_tag("solver"); - inline constexpr auto integrator = make_tag("integrator"); - inline constexpr auto mapping = make_tag("mapping"); - inline constexpr auto utils = make_tag("utils"); - inline constexpr auto mfem_operators = make_tag("operators"); - inline constexpr auto initialization = make_tag("initialization"); - inline constexpr auto accuracy = make_tag("accuracy"); - inline constexpr auto closure = make_tag("closure"); - inline constexpr auto kernels = make_tag("kernels"); - inline constexpr auto surface = make_tag("surface"); - inline constexpr auto model = make_tag("model"); + inline constexpr auto geometry = make_tag("geometry"); + inline constexpr auto physics = make_tag("physics"); + inline constexpr auto unit = make_tag("unit"); + inline constexpr auto mesh = make_tag("mesh"); + inline constexpr auto integration = make_tag("integration"); + inline constexpr auto solver = make_tag("solver"); + inline constexpr auto integrator = make_tag("integrator"); + inline constexpr auto mapping = make_tag("mapping"); + inline constexpr auto utils = make_tag("utils"); + inline constexpr auto mfem_operators = make_tag("operators"); + inline constexpr auto initialization = make_tag("initialization"); + inline constexpr auto accuracy = make_tag("accuracy"); + inline constexpr auto closure = make_tag("closure"); + inline constexpr auto kernels = make_tag("kernels"); + inline constexpr auto surface = make_tag("surface"); + inline constexpr auto model = make_tag("model"); - inline constexpr auto field = sub_tag(mesh & physics, "field"); + inline constexpr auto field = sub_tag(mesh & physics, "field"); - inline constexpr auto legacy_comparison = make_tag("legacy_comparison"); - inline constexpr auto pressure = sub_tag(physics, "pressure"); + inline constexpr auto legacy_comparison = make_tag("legacy_comparison"); + inline constexpr auto pressure = sub_tag(physics, "pressure"); - inline constexpr auto hydro = sub_tag(physics, "hydro"); - inline constexpr auto jacobian = sub_tag(integration & physics, "jacobian"); - inline constexpr auto residuals = sub_tag(integration & physics, "residuals"); - inline constexpr auto volume = sub_tag(mesh & geometry, "volume"); - inline constexpr auto quadrature = sub_tag(mesh & geometry & solver, "quadrature"); - inline constexpr auto convergence = sub_tag(solver, "convergence"); - inline constexpr auto transformations = sub_tag(mesh & geometry, "transformations"); + inline constexpr auto hydro = sub_tag(physics, "hydro"); + inline constexpr auto jacobian = sub_tag(integration & physics, "jacobian"); + inline constexpr auto residuals = sub_tag(integration & physics, "residuals"); + inline constexpr auto volume = sub_tag(mesh & geometry, "volume"); + inline constexpr auto quadrature = sub_tag(mesh & geometry & solver, "quadrature"); + inline constexpr auto convergence = sub_tag(solver, "convergence"); + inline constexpr auto transformations = sub_tag(mesh & geometry, "transformations"); - inline constexpr auto h_refinement = sub_tag(mesh & convergence, "h_refinement"); - inline constexpr auto p_refinement = sub_tag(mesh & convergence, "p_refinement"); + inline constexpr auto h_refinement = sub_tag(mesh & convergence, "h_refinement"); + inline constexpr auto p_refinement = sub_tag(mesh & convergence, "p_refinement"); - inline constexpr auto analytic_comparison = sub_tag(solver & physics & residuals, "analytic_comparison"); - inline constexpr auto self_consistency = sub_tag(solver & physics, "self_consistency"); + inline constexpr auto analytic_comparison = sub_tag(solver & physics & residuals, "analytic_comparison"); + inline constexpr auto self_consistency = sub_tag(solver & physics, "self_consistency"); - inline constexpr auto centrifugal = sub_tag(solver & physics, "centrifugal"); - inline constexpr auto advection = sub_tag(solver & physics, "advection"); - inline constexpr auto coriolis = sub_tag(solver & physics, "coriolis"); - inline constexpr auto gravity = sub_tag(solver & physics, "gravity"); - inline constexpr auto enthalpy = sub_tag(solver & physics, "enthalpy"); - inline constexpr auto barotrope = sub_tag(physics, "barotrope"); - inline constexpr auto mass_continuity = sub_tag(solver & physics, "mass_continuity"); - inline constexpr auto pressure_gradient = sub_tag(solver & physics, "pressure_gradient"); - inline constexpr auto viscosity = sub_tag(solver & physics, "viscosity"); + inline constexpr auto centrifugal = sub_tag(solver & physics, "centrifugal"); + inline constexpr auto advection = sub_tag(solver & physics, "advection"); + inline constexpr auto coriolis = sub_tag(solver & physics, "coriolis"); + inline constexpr auto gravity = sub_tag(solver & physics, "gravity"); + inline constexpr auto enthalpy = sub_tag(solver & physics, "enthalpy"); + inline constexpr auto barotrope = sub_tag(physics, "barotrope"); + inline constexpr auto mass_continuity = sub_tag(solver & physics, "mass_continuity"); + inline constexpr auto pressure_gradient = sub_tag(solver & physics, "pressure_gradient"); + inline constexpr auto viscosity = sub_tag(solver & physics, "viscosity"); - inline constexpr auto compactification = sub_tag(mesh & mapping, "compactification"); - inline constexpr auto kelvin = sub_tag(compactification, "kelvin"); + inline constexpr auto compactification = sub_tag(mesh & mapping, "compactification"); + inline constexpr auto kelvin = sub_tag(compactification, "kelvin"); - inline constexpr auto prepared = sub_tag(solver & physics, "prepared"); - inline constexpr auto contexts = sub_tag(solver, "contexts"); + inline constexpr auto prepared = sub_tag(solver & physics, "prepared"); + inline constexpr auto contexts = sub_tag(solver, "contexts"); - inline constexpr auto domain = sub_tag(mesh, "domain"); + inline constexpr auto domain = sub_tag(mesh, "domain"); + + // Canonical gravity-suite tags. These intentionally compose leaf tags + // exactly once so Catch2 output remains useful and free of repeated + // [solver]/[physics] entries inherited from older composite tags. + inline constexpr auto gravity_unit = gravity & unit; + inline constexpr auto gravity_integration = gravity & integration; + inline constexpr auto gravity_operator = gravity & mfem_operators; + inline constexpr auto gravity_prepared = gravity & make_tag("prepared"); + inline constexpr auto gravity_context = gravity & make_tag("context"); + inline constexpr auto gravity_kernel = gravity & kernels; + inline constexpr auto gravity_accuracy = gravity & accuracy; + inline constexpr auto gravity_legacy = gravity & legacy_comparison; + inline constexpr auto gravity_operator_unit = gravity_operator & unit; + inline constexpr auto gravity_operator_integration = gravity_operator & integration; + inline constexpr auto gravity_operator_convergence = gravity_operator & integration & make_tag("convergence"); + inline constexpr auto gravity_analytic = gravity & integration & make_tag("analytic_comparison"); + inline constexpr auto gravity_consistency = gravity & integration & make_tag("self_consistency"); + inline constexpr auto gravity_prepared_jacobian = gravity_prepared & integration & make_tag("jacobian"); + inline constexpr auto gravity_prepared_unit = gravity_prepared & unit; + inline constexpr auto gravity_prepared_jacobian_accuracy = gravity_prepared_jacobian & accuracy; + inline constexpr auto gravity_kernel_accuracy = gravity_kernel & accuracy; + inline constexpr auto gravity_kernel_integration = gravity_kernel & integration; + inline constexpr auto gravity_kernel_convergence = gravity_kernel & integration & make_tag("convergence"); + inline constexpr auto gravity_analytic_accuracy = gravity_analytic & accuracy; + inline constexpr auto gravity_analytic_initialization = gravity_analytic & initialization; + inline constexpr auto gravity_consistency_initialization = gravity_consistency & initialization; + + inline constexpr auto barotrope_prepared = barotrope & solver & make_tag("prepared"); + inline constexpr auto barotrope_prepared_jacobian = barotrope_prepared & integration & make_tag("jacobian"); + inline constexpr auto barotrope_context = barotrope & solver & make_tag("context"); + inline constexpr auto barotrope_context_integration = barotrope_context & integration; + inline constexpr auto barotrope_prepared_analytic = + barotrope_prepared & integration & make_tag("analytic_comparison"); + inline constexpr auto barotrope_prepared_jacobian_accuracy = barotrope_prepared_jacobian & accuracy; + inline constexpr auto barotrope_prepared_jacobian_geometry = barotrope_prepared_jacobian & geometry; + inline constexpr auto barotrope_prepared_jacobian_unit = barotrope_prepared_jacobian & unit; + + // Canonical hydrostatic-suite tags. The leaf tags are composed directly + // so inherited [physics]/[solver] tags appear only once. + inline constexpr auto barotrope_hydrostatic = barotrope & solver & make_tag("hydro"); + inline constexpr auto barotrope_hydrostatic_context = barotrope_hydrostatic & make_tag("context"); + inline constexpr auto barotrope_hydrostatic_prepared = barotrope_hydrostatic & make_tag("prepared"); + inline constexpr auto barotrope_hydrostatic_prepared_residual = + barotrope_hydrostatic_prepared & integration & make_tag("residual"); + inline constexpr auto barotrope_hydrostatic_prepared_jacobian = + barotrope_hydrostatic_prepared & integration & make_tag("jacobian"); + inline constexpr auto barotrope_hydrostatic_prepared_analytic = + barotrope_hydrostatic_prepared & integration & make_tag("analytic_comparison"); + + inline constexpr auto barotrope_mass_normalization = + barotrope & solver & make_tag("mass_normalization"); + inline constexpr auto barotrope_mass_normalization_context = + barotrope_mass_normalization & make_tag("context"); + inline constexpr auto barotrope_mass_normalization_prepared = + barotrope_mass_normalization & make_tag("prepared"); + inline constexpr auto barotrope_mass_normalization_jacobian = + barotrope_mass_normalization_prepared & integration & make_tag("jacobian"); + inline constexpr auto barotrope_mass_normalization_analytic = + barotrope_mass_normalization_prepared & integration & make_tag("analytic_comparison"); + + inline constexpr auto rotation_prepared = centrifugal & make_tag("prepared"); + inline constexpr auto rotation_context = centrifugal & make_tag("context"); + inline constexpr auto rotation_analytic = centrifugal & integration & make_tag("analytic_comparison"); + inline constexpr auto rotation_context_unit = rotation_context & unit; + inline constexpr auto rotation_prepared_unit = rotation_prepared & unit; + inline constexpr auto rotation_prepared_jacobian = rotation_prepared & integration & make_tag("jacobian"); + inline constexpr auto rotation_prepared_jacobian_accuracy = rotation_prepared_jacobian & accuracy; + inline constexpr auto rotation_kernel_accuracy = centrifugal & kernels & accuracy; + inline constexpr auto rotation_analytic_unit = rotation_analytic & unit; + inline constexpr auto rotation_analytic_accuracy = rotation_analytic & accuracy; + inline constexpr auto rotation_analytic_accuracy_geometry = rotation_analytic_accuracy & geometry; } // namespace tags