457 lines
30 KiB
C++
457 lines
30 KiB
C++
#pragma once
|
|
|
|
// Include after import mean_field. These observations use the already prepared
|
|
// base and private grid functions; they do not prepare a trial or alter a solver
|
|
// buffer. Scalar Eulerian variations are local, broken derivatives for DG
|
|
// fields. They do not contain jump distributions and cannot establish a gauge.
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <cstring>
|
|
#include <filesystem>
|
|
#include <fstream>
|
|
#include <iomanip>
|
|
#include <limits>
|
|
#include <memory>
|
|
#include <stdexcept>
|
|
#include <string>
|
|
#include <string_view>
|
|
#include <utility>
|
|
|
|
#include <mfem.hpp>
|
|
#include <mpi.h>
|
|
|
|
namespace experiment::plateau_diagnostics {
|
|
namespace shape_field_physical_detail {
|
|
using DomainSchema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema;
|
|
|
|
// DomainMapper caches reference geometry even through its const interface.
|
|
// A private instance prevents the observation from warming the solver's cache.
|
|
// Its exterior policy deliberately has no numerical behavior: all observations
|
|
// below must be interior under both this map and the prepared production map.
|
|
class RejectExteriorMap final : public mean_field::mapping::compactification::ExteriorDomainMap {
|
|
public:
|
|
mean_field::mapping::MappingStatus Evaluate(
|
|
const mean_field::mapping::compactification::ExteriorMapInput &,
|
|
mean_field::mapping::compactification::ExteriorMapResult &) const override {
|
|
throw std::logic_error("Shape-field physical diagnostic cannot evaluate the exterior domain.");
|
|
}
|
|
mean_field::mapping::MappingStatus EvaluateVariation(
|
|
const mean_field::mapping::compactification::ExteriorMapInput &,
|
|
const mean_field::mapping::compactification::ExteriorMapResult &,
|
|
const mean_field::mapping::compactification::ExteriorMapDirection &,
|
|
mean_field::mapping::compactification::ExteriorMapVariation &) const override {
|
|
throw std::logic_error("Shape-field physical diagnostic cannot vary the exterior domain.");
|
|
}
|
|
std::string_view GetName() const noexcept override { return "diagnostic_interior_only"; }
|
|
};
|
|
|
|
inline void Require(bool condition, const std::string &message) {
|
|
if (!condition) throw std::runtime_error("Shape-field physical diagnostic: " + message);
|
|
}
|
|
|
|
inline bool Finite(const mfem::Vector &vector) {
|
|
for (int i = 0; i < vector.Size(); ++i)
|
|
if (!std::isfinite(vector(i))) return false;
|
|
return true;
|
|
}
|
|
|
|
inline bool Identical(const mfem::Vector &left, const mfem::Vector &right) {
|
|
return left.Size() == right.Size() &&
|
|
(left.Size() == 0 || std::memcmp(left.GetData(), right.GetData(),
|
|
sizeof(mfem::real_t) * left.Size()) == 0);
|
|
}
|
|
|
|
inline std::string Quote(std::string_view text) {
|
|
std::string result = "\"";
|
|
for (const char ch : text) {
|
|
if (ch == '"') result += '"';
|
|
result += ch;
|
|
}
|
|
return result + '"';
|
|
}
|
|
|
|
inline std::ofstream Output(const std::filesystem::path &directory, std::string_view suffix) {
|
|
const auto path = directory / ("weak_shape_physical_" + std::string(suffix) + ".csv");
|
|
Require(!std::filesystem::exists(path), "output already exists: " + path.string());
|
|
std::ofstream output(path);
|
|
output.exceptions(std::ios::failbit | std::ios::badbit);
|
|
output << std::setprecision(std::numeric_limits<long double>::max_digits10);
|
|
return output;
|
|
}
|
|
|
|
struct Block { int offset; int size; };
|
|
|
|
template <class Blocks>
|
|
Block RequireBlock(const Blocks &blocks, std::string_view name, int size) {
|
|
for (const auto &block : blocks) {
|
|
if (block.stableId != name) continue;
|
|
Require(block.offset >= 0 && block.size >= 0 && block.offset <= size - block.size,
|
|
"invalid root block " + std::string(name));
|
|
return {block.offset, block.size};
|
|
}
|
|
throw std::invalid_argument("Shape-field physical diagnostic requires root block " + std::string(name));
|
|
}
|
|
|
|
template <class Field, class Blocks>
|
|
double Scatter(const Blocks &blocks, std::string_view name, const mfem::Vector &values,
|
|
mfem::ParGridFunction &field) {
|
|
const auto block = RequireBlock(blocks, name, values.Size());
|
|
const auto adapter = mean_field::field::make_field_dof_grid_function_adapter<Field, DomainSchema>(
|
|
*field.ParFESpace());
|
|
Require(block.size == adapter.dof_map().reduced_size(), "field adapter size disagrees with " + std::string(name));
|
|
mfem::Vector reduced(block.size);
|
|
for (int i = 0; i < block.size; ++i) reduced(i) = values(block.offset + i);
|
|
adapter.scatter(reduced, field);
|
|
const auto gathered = adapter.gather(field);
|
|
double maximumError = 0.;
|
|
for (int i = 0; i < block.size; ++i)
|
|
maximumError = std::max(maximumError, std::abs(gathered(i) - reduced(i)));
|
|
Require(Finite(field) && std::isfinite(maximumError), "nonfinite field reconstruction");
|
|
return maximumError;
|
|
}
|
|
|
|
// Every entry is a raw additive physical-volume integral. Keeping the signed
|
|
// cross term lets an independent reader verify E^2 = D^2 + T^2 - 2 D.T.
|
|
struct ScalarIntegrals {
|
|
long double base{}, pullback{}, transport{}, eulerian{};
|
|
long double baseSquared{}, pullbackSquared{}, transportSquared{}, eulerianSquared{}, cross{};
|
|
|
|
void Add(double weight, double f, double delta, double advection) {
|
|
const long double w = weight, b = f, d = delta, t = advection, e = d - t;
|
|
base += w * b;
|
|
pullback += w * d;
|
|
transport += w * t;
|
|
eulerian += w * e;
|
|
baseSquared += w * b * b;
|
|
pullbackSquared += w * d * d;
|
|
transportSquared += w * t * t;
|
|
eulerianSquared += w * e * e;
|
|
cross += w * d * t;
|
|
}
|
|
|
|
ScalarIntegrals &operator+=(const ScalarIntegrals &other) {
|
|
base += other.base; pullback += other.pullback; transport += other.transport; eulerian += other.eulerian;
|
|
baseSquared += other.baseSquared; pullbackSquared += other.pullbackSquared;
|
|
transportSquared += other.transportSquared; eulerianSquared += other.eulerianSquared; cross += other.cross;
|
|
return *this;
|
|
}
|
|
|
|
void Write(std::ostream &output) const {
|
|
for (const auto value : {base, pullback, transport, eulerian, baseSquared, pullbackSquared,
|
|
transportSquared, eulerianSquared, cross})
|
|
Require(std::isfinite(value), "nonfinite integrated scalar quantity");
|
|
output << ',' << base << ',' << pullback << ',' << transport << ',' << eulerian << ',' << baseSquared
|
|
<< ',' << pullbackSquared << ',' << transportSquared << ',' << eulerianSquared << ',' << cross;
|
|
}
|
|
};
|
|
|
|
inline constexpr std::string_view scalarColumns =
|
|
",base_integral,pullback_integral,transport_integral,eulerian_integral,base_squared_integral,"
|
|
"pullback_squared_integral,transport_squared_integral,eulerian_squared_integral,pullback_transport_cross_integral";
|
|
} // namespace shape_field_physical_detail
|
|
|
|
template <class State>
|
|
void MeasureShapeFieldPhysical(State &state, const mean_field::fem::FEM &fem,
|
|
const mfem::Vector &physicalBase, const mfem::Vector &physicalDirection,
|
|
const std::filesystem::path &output) {
|
|
namespace d = shape_field_physical_detail;
|
|
using namespace mean_field;
|
|
const auto &problem = state.Problem();
|
|
int ranks = 0;
|
|
MPI_Comm_size(problem.GetCommunicator(), &ranks);
|
|
d::Require(ranks == 1, "exactly one MPI rank is required");
|
|
d::Require(fem.mesh && fem.displacementFes && fem.surfaceDeformationFes && fem.densityFes &&
|
|
fem.enthalpyFes && fem.gravityPotentialFes && fem.gravityFluxFes && fem.compactificationCoordinate,
|
|
"incomplete finite-element spaces");
|
|
d::Require(fem.mesh->SpaceDimension() == 3 && fem.mesh->Dimension() == 3,
|
|
"three-dimensional volume elements are required");
|
|
d::Require(state.normalizedOperator->IsPrepared() && problem.GetPhysicalOperator().IsPrepared(),
|
|
"the base operator is not prepared");
|
|
d::Require(physicalBase.Size() == problem.StateSize() && physicalDirection.Size() == physicalBase.Size() &&
|
|
d::Finite(physicalBase) && d::Finite(physicalDirection) &&
|
|
d::Identical(physicalBase, state.normalizedOperator->GetPhysicalState()),
|
|
"the finite supplied base must exactly equal the prepared physical state");
|
|
|
|
// Capture what a caller could otherwise mistake for a silent preparation.
|
|
// The only production action below is the immutable deformation Jacobian.
|
|
const auto generation = problem.GetPreparationGeneration();
|
|
const mfem::Vector accepted(state.acceptedNormalizedState), acceptedResidual(state.acceptedNormalizedResidual),
|
|
acceptedPhysical(state.AcceptedPhysicalState()), trial(state.trialNormalizedState),
|
|
trialResidual(state.trialNormalizedResidual), candidatePhysical(state.candidatePhysicalState),
|
|
preparedResidual(state.normalizedOperator->GetPhysicalResidual()),
|
|
stateFactors(state.normalizedOperator->GetNormalization().StateFactors()),
|
|
residualFactors(state.normalizedOperator->GetNormalization().ResidualFactors());
|
|
const auto &physical = problem.GetPhysicalOperator();
|
|
const mfem::Vector baseDisplacement(physical.GetGeneratedVolumeDisplacement());
|
|
mfem::Vector volumeDirection(fem.displacementFes->GetTrueVSize());
|
|
problem.BuildVolumeDisplacementDirection(physicalDirection, volumeDirection);
|
|
d::Require(volumeDirection.Size() == fem.displacementFes->GetTrueVSize() && d::Finite(volumeDirection),
|
|
"invalid generated volume direction");
|
|
|
|
mfem::ParGridFunction rho(fem.densityFes.get()), drho(fem.densityFes.get()),
|
|
h(fem.enthalpyFes.get()), dh(fem.enthalpyFes.get()),
|
|
phi(fem.gravityPotentialFes.get()), dphi(fem.gravityPotentialFes.get()),
|
|
displacement(fem.displacementFes.get()), velocity(fem.displacementFes.get()),
|
|
affine(fem.displacementFes.get());
|
|
const auto blocks = problem.GetManifest().valueBlocks();
|
|
const std::array<double, 6> scatterErrors{
|
|
d::Scatter<field::Density>(blocks, "density", physicalBase, rho),
|
|
d::Scatter<field::Density>(blocks, "density", physicalDirection, drho),
|
|
d::Scatter<field::Enthalpy>(blocks, "specific_enthalpy", physicalBase, h),
|
|
d::Scatter<field::Enthalpy>(blocks, "specific_enthalpy", physicalDirection, dh),
|
|
d::Scatter<field::Gravity>(blocks, "gravity_potential", physicalBase, phi),
|
|
d::Scatter<field::Gravity>(blocks, "gravity_potential", physicalDirection, dphi)};
|
|
displacement.SetFromTrueDofs(baseDisplacement);
|
|
velocity.SetFromTrueDofs(volumeDirection);
|
|
mfem::VectorFunctionCoefficient identity(3, [](const mfem::Vector &x, mfem::Vector &value) { value = x; });
|
|
affine.ProjectCoefficient(identity);
|
|
d::Require(d::Finite(displacement) && d::Finite(velocity) && d::Finite(affine), "nonfinite private displacement field");
|
|
const auto &productionMapper = problem.GetDiscretization().domainMapper();
|
|
const auto productionCacheBefore = productionMapper.GetGeometryCacheStatistics();
|
|
mapping::DomainMapper privateMapper(
|
|
utils::DomainMapperOptions{.dimension = 3,
|
|
.vacuum_element_attribute = d::DomainSchema::template material_attribute<utils::domain::Vacuum>()},
|
|
std::make_unique<d::RejectExteriorMap>());
|
|
privateMapper.SetGeometryCacheEnabled(false);
|
|
mapping::GridFunctionMappingEvaluator mappingEvaluator(privateMapper, displacement, *fem.compactificationCoordinate);
|
|
mappingEvaluator.InvalidateCache();
|
|
|
|
// The default surface descriptor erases its center. Checking both actual
|
|
// displacement traces prevents a different radial prescription from being
|
|
// silently interpreted in the zero-center basis used for this CSV.
|
|
const auto &domain = physical.GetDomainDeformation();
|
|
mfem::Vector center(3); center = 0.;
|
|
const deformation::NodalRadialSurface prescription(std::move(center));
|
|
d::Require(domain.descriptor().surfaceDeformation == prescription.descriptor() &&
|
|
fem.surfaceDeformationFes->GetVDim() == 1 && fem.displacementFes->GetVDim() == 3 &&
|
|
fem.displacementFes->GetOrdering() == mfem::Ordering::byNODES &&
|
|
fem.displacementFes->FEColl() == fem.surfaceDeformationFes->FEColl() &&
|
|
fem.displacementFes->GetTrueVSize() == 3 * fem.surfaceDeformationFes->GetTrueVSize(),
|
|
"the default nodal radial surface and its vector companion are required");
|
|
const auto radial = deformation::compileSurfaceDeformationPrescription(prescription,
|
|
deformation::SurfaceDeformationCompilationContext{
|
|
*fem.surfaceDeformationFes, field::make_stellar_surface_scalar_dof_map(*fem.surfaceDeformationFes)});
|
|
const auto shape = d::RequireBlock(blocks, "surface_deformation", physicalBase.Size());
|
|
d::Require(radial.parameterCount() == shape.size && shape.size == domain.parameterCount(),
|
|
"surface parameter counts disagree");
|
|
const int scalarSize = fem.surfaceDeformationFes->GetTrueVSize();
|
|
double traceError = 0., traceScale = 1.;
|
|
for (int i = 0; i < shape.size; ++i)
|
|
for (int component = 0; component < 3; ++component) {
|
|
const int dof = radial.surfaceDofMap().volume_true_dof(i) + component * scalarSize;
|
|
const double direction = radial.radialDirection(i, component);
|
|
const double expectedBase = physicalBase(shape.offset + i) * direction;
|
|
const double expectedDirection = physicalDirection(shape.offset + i) * direction;
|
|
traceError = std::max({traceError, std::abs(baseDisplacement(dof) - expectedBase),
|
|
std::abs(volumeDirection(dof) - expectedDirection)});
|
|
traceScale = std::max({traceScale, std::abs(expectedBase), std::abs(expectedDirection)});
|
|
}
|
|
d::Require(std::isfinite(traceError) && traceError <= 512. * std::numeric_limits<double>::epsilon() * traceScale,
|
|
"actual generated surface trace disagrees with the reconstructed radial basis");
|
|
|
|
std::filesystem::create_directories(output);
|
|
auto metadata = d::Output(output, "metadata");
|
|
auto elements = d::Output(output, "elements");
|
|
auto geometry = d::Output(output, "geometry");
|
|
auto affineOutput = d::Output(output, "affine");
|
|
auto surface = d::Output(output, "surface_nodes");
|
|
auto summary = d::Output(output, "summary");
|
|
metadata << "key,value\n";
|
|
const auto meta = [&](std::string_view key, std::string_view value) {
|
|
metadata << d::Quote(key) << ',' << d::Quote(value) << '\n';
|
|
};
|
|
meta("schema", "weak_shape_physical_v1");
|
|
meta("scope", "stellar interior only; prepared base; physical direction supplied; no trial preparation");
|
|
meta("scalar_formula", "delta_f_E = delta_f_at_fixed_reference - grad_physical_f dot volume_displacement_direction");
|
|
meta("dg_derivatives", "density and potential gradients are element-local broken derivatives; jumps excluded");
|
|
meta("interpretation", "local shape-field compensation diagnostic; no nullspace or pure-gauge certification");
|
|
meta("quadrature", "fixed integration-rule orders 10 and 14, independent of field order; compare for integration sensitivity");
|
|
meta("radius", "Euclidean distance from zero; per-element min/max are quadrature-sample extrema");
|
|
meta("surface_samples", "nodal values without area weights; not harmonic coefficients or surface integrals");
|
|
meta("affine_test", "projection of reference physical coordinate x into displacement space; independent of the actual domain lift");
|
|
meta("affine_weights", "reference physical mesh volume; gradient is with respect to reference physical coordinates");
|
|
meta("physical_weights", "integration_point.weight * element_transformation.Weight * mapping_determinant");
|
|
meta("mapping", "private uncached DomainMapper; same interior x+u implementation; exterior evaluation forbidden");
|
|
meta("scalar_integrals", "raw additive integrals; transport = grad_physical_f dot v; no distributional interface contribution");
|
|
meta("gravity_vector", "not evaluated; RT Piola variation is outside this scalar diagnostic");
|
|
meta("reference_mesh_geometry_collection", fem.mesh->GetNodes()
|
|
? fem.mesh->GetNodes()->FESpace()->FEColl()->Name() : "linear_mesh_without_nodal_grid_function");
|
|
const std::array<std::pair<std::string_view, const mfem::ParFiniteElementSpace *>, 6> spaces{{
|
|
{"density", fem.densityFes.get()}, {"specific_enthalpy", fem.enthalpyFes.get()},
|
|
{"gravity_potential", fem.gravityPotentialFes.get()}, {"displacement", fem.displacementFes.get()},
|
|
{"surface_deformation", fem.surfaceDeformationFes.get()}, {"gravity_gradient", fem.gravityFluxFes.get()}}};
|
|
for (const auto &[name, space] : spaces) {
|
|
meta(std::string(name) + "_collection", space->FEColl()->Name());
|
|
meta(std::string(name) + "_true_size", std::to_string(space->GetTrueVSize()));
|
|
}
|
|
meta("density_family_order", std::to_string(field::Density::Scalar::familyOrder));
|
|
meta("enthalpy_family_order", std::to_string(field::Enthalpy::Scalar::familyOrder));
|
|
meta("potential_family_order", std::to_string(field::Gravity::Potential::familyOrder));
|
|
meta("gravity_flux_family_order", std::to_string(field::Gravity::Flux::familyOrder));
|
|
meta("displacement_family_order", std::to_string(field::Displacement::Vector::familyOrder));
|
|
|
|
elements << "quadrature,element,attribute,geometry,field,gradient_kind,points,physical_volume"
|
|
<< d::scalarColumns << '\n';
|
|
geometry << "quadrature,element,attribute,geometry,points,reference_mesh_geometry_order,density_element_order,enthalpy_element_order,"
|
|
"potential_element_order,displacement_element_order,reference_radius_min,reference_radius_max,"
|
|
"physical_radius_min,physical_radius_max,physical_volume,base_displacement_squared_integral,"
|
|
"direction_displacement_squared_integral\n";
|
|
affineOutput << "quadrature,element,points,reference_volume,value_error_squared_integral,"
|
|
"reference_value_squared_integral,gradient_error_squared_integral,reference_gradient_squared_integral\n";
|
|
surface << "parameter,volume_scalar_true_dof,unit_radial_x,unit_radial_y,unit_radial_z,reference_radius,"
|
|
"physical_radius,base_parameter,physical_direction_parameter\n";
|
|
for (int i = 0; i < shape.size; ++i) {
|
|
const double radius = radial.referenceRadius(i), base = physicalBase(shape.offset + i);
|
|
d::Require(std::isfinite(radius) && radius > 0. && radius + base > 0., "invalid radial surface radius");
|
|
surface << i << ',' << radial.surfaceDofMap().volume_true_dof(i);
|
|
for (int component = 0; component < 3; ++component) surface << ',' << radial.radialDirection(i, component);
|
|
surface << ',' << radius << ',' << radius + base << ',' << base << ',' << physicalDirection(shape.offset + i) << '\n';
|
|
}
|
|
summary << "quadrature,field,metric,value\n";
|
|
const auto measure = [&](int q, std::string_view fieldName, std::string_view metric, long double value) {
|
|
d::Require(std::isfinite(value), "nonfinite summary " + std::string(metric));
|
|
summary << q << ',' << fieldName << ',' << metric << ',' << value << '\n';
|
|
};
|
|
const std::array<std::string_view, 3> names{"density", "specific_enthalpy", "gravity_potential"};
|
|
const std::array<mfem::ParGridFunction *, 3> fields{&rho, &h, &phi}, directions{&drho, &dh, &dphi};
|
|
mfem::Vector referenceGradient(3), physicalGradient(3), u(3), v(3), x(3), affineValue(3);
|
|
mfem::DenseMatrix affineGradient(3);
|
|
mapping::MappingPointContext mapped;
|
|
for (const int quadrature : {10, 14}) {
|
|
std::array<d::ScalarIntegrals, 3> totals;
|
|
long double totalVolume = 0., totalU2 = 0., totalV2 = 0., affineVolume = 0., affineError2 = 0.,
|
|
affineReference2 = 0., affineGradientError2 = 0.;
|
|
int stellarElements = 0;
|
|
for (int element = 0; element < fem.mesh->GetNE(); ++element) {
|
|
const int attribute = fem.mesh->GetAttribute(element);
|
|
if (!d::DomainSchema::template attribute_belongs_to<utils::domain::Stellar>(attribute)) continue;
|
|
++stellarElements;
|
|
auto *transformation = fem.mesh->GetElementTransformation(element);
|
|
d::Require(!productionMapper.IsCompactifiedElement(*transformation),
|
|
"production mapper classifies a stellar element as compactified");
|
|
const auto &rule = mfem::IntRules.Get(transformation->GetGeometryType(), quadrature);
|
|
std::array<d::ScalarIntegrals, 3> integrals;
|
|
long double volume = 0., u2 = 0., v2 = 0., referenceVolume = 0., affineE2 = 0., affineR2 = 0., affineGE2 = 0.;
|
|
double rmin = std::numeric_limits<double>::infinity(), rmax = 0.,
|
|
physicalRmin = std::numeric_limits<double>::infinity(), physicalRmax = 0.;
|
|
for (int point = 0; point < rule.GetNPoints(); ++point) {
|
|
const auto &ip = rule.IntPoint(point);
|
|
d::Require(mappingEvaluator.EvaluatePoint(*transformation, ip, mapped) == mapping::MappingStatus::valid &&
|
|
!mapped.compactified && std::isfinite(mapped.mapping_determinant) && mapped.mapping_determinant > 0.,
|
|
"invalid or compactified stellar mapping");
|
|
transformation->SetIntPoint(&ip);
|
|
const double referenceWeight = ip.weight * transformation->Weight();
|
|
const double weight = referenceWeight * mapped.mapping_determinant;
|
|
d::Require(std::isfinite(referenceWeight) && referenceWeight > 0. && std::isfinite(weight) && weight > 0.,
|
|
"invalid quadrature weight");
|
|
displacement.GetVectorValue(element, ip, u);
|
|
velocity.GetVectorValue(element, ip, v);
|
|
d::Require(d::Finite(u) && d::Finite(v) && d::Finite(mapped.reference_position) && d::Finite(mapped.physical_position),
|
|
"nonfinite sampled displacement or position");
|
|
const double r = mapped.reference_position.Norml2(), physicalR = mapped.physical_position.Norml2();
|
|
d::Require(std::isfinite(r) && std::isfinite(physicalR), "overflowing sampled radius");
|
|
rmin = std::min(rmin, r); rmax = std::max(rmax, r);
|
|
physicalRmin = std::min(physicalRmin, physicalR); physicalRmax = std::max(physicalRmax, physicalR);
|
|
volume += weight;
|
|
for (int component = 0; component < 3; ++component) {
|
|
u2 += static_cast<long double>(weight) * u(component) * u(component);
|
|
v2 += static_cast<long double>(weight) * v(component) * v(component);
|
|
}
|
|
for (std::size_t fieldIndex = 0; fieldIndex < fields.size(); ++fieldIndex) {
|
|
const double base = fields[fieldIndex]->GetValue(element, ip);
|
|
const double delta = directions[fieldIndex]->GetValue(element, ip);
|
|
transformation->SetIntPoint(&ip);
|
|
fields[fieldIndex]->GetGradient(*transformation, referenceGradient);
|
|
mapped.inverse_mapping_jacobian.MultTranspose(referenceGradient, physicalGradient);
|
|
const double transport = physicalGradient * v;
|
|
d::Require(std::isfinite(base) && std::isfinite(delta) && std::isfinite(transport) && d::Finite(physicalGradient),
|
|
"nonfinite scalar or broken gradient");
|
|
integrals[fieldIndex].Add(weight, base, delta, transport);
|
|
}
|
|
transformation->SetIntPoint(&ip);
|
|
transformation->Transform(ip, x);
|
|
affine.GetVectorValue(element, ip, affineValue);
|
|
affine.GetVectorGradient(*transformation, affineGradient);
|
|
d::Require(d::Finite(x) && d::Finite(affineValue), "nonfinite affine projection sample");
|
|
referenceVolume += referenceWeight;
|
|
for (int row = 0; row < 3; ++row) {
|
|
const long double error = static_cast<long double>(affineValue(row)) - x(row);
|
|
affineE2 += referenceWeight * error * error;
|
|
affineR2 += static_cast<long double>(referenceWeight) * x(row) * x(row);
|
|
for (int column = 0; column < 3; ++column) {
|
|
const long double gradientError = affineGradient(row, column) - (row == column ? 1. : 0.);
|
|
d::Require(std::isfinite(gradientError), "nonfinite affine projection gradient");
|
|
affineGE2 += referenceWeight * gradientError * gradientError;
|
|
}
|
|
}
|
|
}
|
|
for (std::size_t fieldIndex = 0; fieldIndex < fields.size(); ++fieldIndex) {
|
|
elements << quadrature << ',' << element << ',' << attribute << ',' << int(transformation->GetGeometryType())
|
|
<< ',' << names[fieldIndex] << ',' << (fieldIndex == 1 ? "continuous_H1" : "broken_DG")
|
|
<< ',' << rule.GetNPoints() << ',' << volume;
|
|
integrals[fieldIndex].Write(elements);
|
|
elements << '\n';
|
|
totals[fieldIndex] += integrals[fieldIndex];
|
|
}
|
|
geometry << quadrature << ',' << element << ',' << attribute << ',' << int(transformation->GetGeometryType())
|
|
<< ',' << rule.GetNPoints() << ',' << (fem.mesh->GetNodes()
|
|
? fem.mesh->GetNodes()->FESpace()->GetElementOrder(element) : 1)
|
|
<< ',' << fem.densityFes->GetElementOrder(element)
|
|
<< ',' << fem.enthalpyFes->GetElementOrder(element) << ',' << fem.gravityPotentialFes->GetElementOrder(element)
|
|
<< ',' << fem.displacementFes->GetElementOrder(element) << ',' << rmin << ',' << rmax << ',' << physicalRmin
|
|
<< ',' << physicalRmax << ',' << volume << ',' << u2 << ',' << v2 << '\n';
|
|
affineOutput << quadrature << ',' << element << ',' << rule.GetNPoints() << ',' << referenceVolume << ','
|
|
<< affineE2 << ',' << affineR2 << ',' << affineGE2 << ',' << 3. * referenceVolume << '\n';
|
|
totalVolume += volume; totalU2 += u2; totalV2 += v2;
|
|
affineVolume += referenceVolume; affineError2 += affineE2; affineReference2 += affineR2; affineGradientError2 += affineGE2;
|
|
}
|
|
d::Require(stellarElements > 0 && totalVolume > 0. && affineVolume > 0. && affineReference2 > 0., "empty stellar integrals");
|
|
for (std::size_t i = 0; i < totals.size(); ++i) {
|
|
const auto &integral = totals[i];
|
|
measure(quadrature, names[i], "base_squared_integral", integral.baseSquared);
|
|
measure(quadrature, names[i], "pullback_squared_integral", integral.pullbackSquared);
|
|
measure(quadrature, names[i], "transport_squared_integral", integral.transportSquared);
|
|
measure(quadrature, names[i], "eulerian_squared_integral", integral.eulerianSquared);
|
|
measure(quadrature, names[i], "pullback_transport_cross_integral", integral.cross);
|
|
measure(quadrature, names[i], "pullback_l2_norm", std::sqrt(integral.pullbackSquared));
|
|
measure(quadrature, names[i], "transport_l2_norm", std::sqrt(integral.transportSquared));
|
|
measure(quadrature, names[i], "eulerian_l2_norm", std::sqrt(integral.eulerianSquared));
|
|
measure(quadrature, names[i], "squared_norm_identity_error",
|
|
integral.eulerianSquared - integral.pullbackSquared - integral.transportSquared + 2. * integral.cross);
|
|
}
|
|
measure(quadrature, "geometry", "stellar_elements", stellarElements);
|
|
measure(quadrature, "geometry", "physical_volume", totalVolume);
|
|
measure(quadrature, "geometry", "base_displacement_squared_integral", totalU2);
|
|
measure(quadrature, "geometry", "direction_displacement_squared_integral", totalV2);
|
|
measure(quadrature, "affine", "reference_volume", affineVolume);
|
|
measure(quadrature, "affine", "value_error_squared_integral", affineError2);
|
|
measure(quadrature, "affine", "reference_value_squared_integral", affineReference2);
|
|
measure(quadrature, "affine", "gradient_error_squared_integral", affineGradientError2);
|
|
measure(quadrature, "affine", "reference_gradient_squared_integral", 3. * affineVolume);
|
|
measure(quadrature, "affine", "relative_value_l2_error", std::sqrt(affineError2 / affineReference2));
|
|
measure(quadrature, "affine", "relative_gradient_l2_error", std::sqrt(affineGradientError2 / (3. * affineVolume)));
|
|
}
|
|
measure(0, "reconstruction", "maximum_scatter_roundtrip_error", *std::max_element(scatterErrors.begin(), scatterErrors.end()));
|
|
measure(0, "reconstruction", "maximum_surface_trace_error", traceError);
|
|
const auto productionCacheAfter = productionMapper.GetGeometryCacheStatistics();
|
|
const bool cacheUnchanged = productionCacheBefore.hits == productionCacheAfter.hits &&
|
|
productionCacheBefore.misses == productionCacheAfter.misses && productionCacheBefore.points == productionCacheAfter.points;
|
|
const bool unchanged = cacheUnchanged &&
|
|
generation == problem.GetPreparationGeneration() && state.normalizedOperator->IsPrepared() &&
|
|
d::Identical(physicalBase, state.normalizedOperator->GetPhysicalState()) &&
|
|
d::Identical(preparedResidual, state.normalizedOperator->GetPhysicalResidual()) &&
|
|
d::Identical(accepted, state.acceptedNormalizedState) && d::Identical(acceptedResidual, state.acceptedNormalizedResidual) &&
|
|
d::Identical(acceptedPhysical, state.AcceptedPhysicalState()) && d::Identical(trial, state.trialNormalizedState) &&
|
|
d::Identical(trialResidual, state.trialNormalizedResidual) && d::Identical(candidatePhysical, state.candidatePhysicalState) &&
|
|
d::Identical(baseDisplacement, physical.GetGeneratedVolumeDisplacement()) &&
|
|
d::Identical(stateFactors, state.normalizedOperator->GetNormalization().StateFactors()) &&
|
|
d::Identical(residualFactors, state.normalizedOperator->GetNormalization().ResidualFactors());
|
|
measure(0, "invariant", "accepted_trial_prepared_normalization_unchanged", unchanged);
|
|
measure(0, "invariant", "production_mapping_cache_unchanged", cacheUnchanged);
|
|
measure(0, "invariant", "trial_preparations", 0);
|
|
measure(0, "invariant", "volume_direction_actions", 1);
|
|
d::Require(unchanged, "solver state or normalization changed during physical observation");
|
|
}
|
|
} // namespace experiment::plateau_diagnostics
|