diff --git a/data/meshes/dfg_cylinder_2D__ref_0.exo b/data/meshes/dfg_cylinder_2D__ref_0.exo new file mode 100644 index 000000000..0cb6d6b7c Binary files /dev/null and b/data/meshes/dfg_cylinder_2D__ref_0.exo differ diff --git a/data/meshes/dfg_cylinder_2D__ref_1.exo b/data/meshes/dfg_cylinder_2D__ref_1.exo new file mode 100644 index 000000000..87ad024c6 Binary files /dev/null and b/data/meshes/dfg_cylinder_2D__ref_1.exo differ diff --git a/data/meshes/dfg_cylinder_2D__ref_2.exo b/data/meshes/dfg_cylinder_2D__ref_2.exo new file mode 100644 index 000000000..02665229a Binary files /dev/null and b/data/meshes/dfg_cylinder_2D__ref_2.exo differ diff --git a/data/meshes/dfg_cylinder_2D__ref_3.exo b/data/meshes/dfg_cylinder_2D__ref_3.exo new file mode 100644 index 000000000..0f9c223e2 Binary files /dev/null and b/data/meshes/dfg_cylinder_2D__ref_3.exo differ diff --git a/src/serac/numerics/CMakeLists.txt b/src/serac/numerics/CMakeLists.txt index 4d7edac16..d2c6ca4bd 100644 --- a/src/serac/numerics/CMakeLists.txt +++ b/src/serac/numerics/CMakeLists.txt @@ -12,6 +12,8 @@ set(numerics_headers solver_config.hpp stdfunction_operator.hpp petsc_solvers.hpp + extrapolation_operator.hpp + fixed_step_bdf_operator.hpp trust_region_solver.hpp dense_petsc.hpp ) diff --git a/src/serac/numerics/equation_solver.cpp b/src/serac/numerics/equation_solver.cpp index 0af69dcff..73cf52045 100644 --- a/src/serac/numerics/equation_solver.cpp +++ b/src/serac/numerics/equation_solver.cpp @@ -89,7 +89,9 @@ class NewtonSolver : public mfem::NewtonSolver { norm = initial_norm = evaluateNorm(x, r); if (print_options.first_and_last && !print_options.iterations) { - mfem::out << "Newton iteration " << std::setw(3) << 0 << " : ||r|| = " << std::setw(13) << norm << "...\n"; + mfem::out << "Newton iteration " << std::setw(3) << 0 + << " : ||r|| = " << std::scientific << std::setw(12) + << std::setprecision(6) << norm << "...\n"; } norm_goal = std::max(rel_tol * initial_norm, abs_tol); @@ -99,9 +101,13 @@ class NewtonSolver : public mfem::NewtonSolver { for (; true; it++) { MFEM_ASSERT(mfem::IsFinite(norm), "norm = " << norm); if (print_options.iterations) { - mfem::out << "Newton iteration " << std::setw(3) << it << " : ||r|| = " << std::setw(13) << norm; + mfem::out << "Newton iteration " << std::setw(3) << it + << " : ||r|| = " << std::scientific << std::setw(12) + << std::setprecision(6) << norm; if (it > 0) { - mfem::out << ", ||r||/||r_0|| = " << std::setw(13) << (initial_norm != 0.0 ? norm / initial_norm : norm); + mfem::out << ", ||r||/||r_0|| = " + << std::scientific << std::setw(12) << std::setprecision(6) + << (initial_norm != 0.0 ? norm / initial_norm : norm); } mfem::out << '\n'; } @@ -193,7 +199,9 @@ class NewtonSolver : public mfem::NewtonSolver { final_norm = norm; if (print_options.summary || (!converged && print_options.warnings) || print_options.first_and_last) { - mfem::out << "Newton: Number of iterations: " << final_iter << '\n' << " ||r|| = " << final_norm << '\n'; + mfem::out << "Newton: Number of iterations: " << final_iter << '\n' + << " ||r|| = " << std::scientific << std::setw(12) + << std::setprecision(6) << final_norm << '\n'; } if (!converged && (print_options.summary || print_options.warnings)) { mfem::out << "Newton: No convergence!\n"; @@ -629,7 +637,9 @@ class TrustRegion : public mfem::NewtonSolver { norm = initial_norm = computeResidual(X, r); norm_goal = std::max(rel_tol * initial_norm, abs_tol); if (print_options.first_and_last && !print_options.iterations) { - mfem::out << "Newton iteration " << std::setw(3) << 0 << " : ||r|| = " << std::setw(13) << norm << "...\n"; + mfem::out << "Newton iteration " << std::setw(3) << 0 + << " : ||r|| = " << std::scientific << std::setw(12) + << std::setprecision(6) << norm << "...\n"; } prec->iterative_mode = false; tr_precond.iterative_mode = false; @@ -659,10 +669,17 @@ class TrustRegion : public mfem::NewtonSolver { for (; true; it++) { MFEM_ASSERT(mfem::IsFinite(norm), "norm = " << norm); if (print_options.iterations) { - mfem::out << "Newton iteration " << std::setw(3) << it << " : ||r|| = " << std::setw(13) << norm; + mfem::out << "Newton iteration " << std::setw(3) << it + << " : ||r|| = " << std::scientific << std::setw(12) + << std::setprecision(6) << norm; if (it > 0) { - mfem::out << ", ||r||/||r_0|| = " << std::setw(13) << (initial_norm != 0.0 ? norm / initial_norm : norm); - mfem::out << ", x_incr = " << std::setw(13) << trResults.d.Norml2(); + mfem::out << ", ||r||/||r_0|| = " << std::scientific << std::setw(12) + << std::setprecision(6) + << (initial_norm != 0.0 ? norm / initial_norm : norm); + + mfem::out << ", x_incr = " + << std::scientific << std::setw(12) << std::setprecision(6) + << trResults.d.Norml2(); } else { mfem::out << ", norm goal = " << std::setw(13) << norm_goal; } @@ -849,8 +866,9 @@ class TrustRegion : public mfem::NewtonSolver { final_norm = norm; if (print_options.summary || (!converged && print_options.warnings) || print_options.first_and_last) { - mfem::out << "Newton: Number of iterations: " << final_iter << '\n' << " ||r|| = " << final_norm << '\n'; - } + mfem::out << "Newton: Number of iterations: " << final_iter << '\n' + << " ||r|| = " << std::scientific << std::setw(12) + << std::setprecision(6) << final_norm << '\n'; } if (!converged && (print_options.summary || print_options.warnings)) { mfem::out << "Newton: No convergence!\n"; } diff --git a/src/serac/numerics/extrapolation_operator.hpp b/src/serac/numerics/extrapolation_operator.hpp new file mode 100644 index 000000000..ed8d4ebad --- /dev/null +++ b/src/serac/numerics/extrapolation_operator.hpp @@ -0,0 +1,176 @@ +// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and +// other Serac Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/** + * @file extrapolation_operator.hpp + * + * @brief Implements an extrapolation operator for time-dependent fields using + * zero-, first-, or second-order extrapolation. + */ + +#pragma once + +#include "serac/infrastructure/logger.hpp" + +namespace serac { + +/** + * @brief Implements a cycle-aware extrapolator for time-dependent fields. + * + * This operator handles extrapolation of any time-dependent field (e.g., pressure, velocity) + * using a compile-time specified maximum order and a runtime cycle count to determine + * how many states are available. + * + * @tparam r_order The desired extrapolation order (0 = zero, 1 = constant, 2 = linear). + */ +template +struct ExtrapolationOperator { + + static_assert( + r_order >= 0 && r_order <= 2, + "Supported extrapolation orders are 0 (zero), 1 (constant), and 2 (linear)." + ); + + /** + * @brief Constructs an extrapolation operator over the given finite element space. + * + * This operator enables time-extrapolation of fields defined in the specified finite element space, + * such as velocity or pressure, using a compile-time extrapolation order. + * + * The extrapolated field φ^{*,k+1} will live in the same space as the provided input states φᵏ, φᵏ⁻¹, etc. + * + * @param space The finite element space in which the extrapolated fields reside. + */ + ExtrapolationOperator( + mfem::ParFiniteElementSpace *space + ) : space_(space) {} + + /** + * @brief Extrapolates a field value φ^{*,k+1} from previous states. + * + * The extrapolated value is computed as: + * φ^{*,r,k+1} = ∑ₘ γₘ · φ^{k−m}, for m = 0..r−1 + * where r = min(r_order, k+1), and γₘ are extrapolation weights. + * + * The input vector previous_states is ordered from newest (φᵏ) to oldest (φ^{k−(r−1)}). + * + * @param previous_states Vector of past states φᵏ, φᵏ⁻¹, … ordered newest to oldest. + * @param cycle Current simulation cycle (corresponds to time step k+1). + * @return Extrapolated field φ^{*,k+1}. + */ + ::serac::FiniteElementState + operator()( + ::std::vector<::serac::FiniteElementState> const &previous_states, + int cycle + ) const { + + SLIC_ERROR_ROOT_IF(cycle <= 0, "Extrapolation requires at least one completed cycle."); + + int const available_order = std::min(cycle, r_order); // Effective order. + + SLIC_ERROR_ROOT_IF( + static_cast(previous_states.size()) != available_order, + axom::fmt::format("Expected at least {} previous states, but got {}.", available_order, previous_states.size()) + ); + + ::serac::FiniteElementState result(*space_); // Initialize from most recent state (φ^k). + + if (available_order == 0) { + result = 0.0; // Zero extrapolation. + } else { + auto gammas = this->compute_gammas(available_order); + result = previous_states[0]; + result *= gammas[0]; // result = γ₀ · φ^k + for (int j = 1; j < available_order; ++j) { + result.Add(gammas[j], previous_states[j]); // result += γⱼ · φ^{k−j} + } + } + + return result; + + } + + /** + * @brief Applies the derivative ∂φ^{*,j} / ∂φᵗ to a vector. + * + * Given a target extrapolated field φ^{*,j} and a source field φᵗ, this method computes: + * ∂φ^{*,j} / ∂φᵗ · v = γₘ · v, + * where m = j − 1 − t, and γₘ is the corresponding extrapolation weight. + * + * If m is outside the valid range [0, r−1], where r = min(r_order, j), the result is zero. + * + * @tparam VectorType A vector-like type supporting assignment and scalar multiplication. + * @param j The target cycle index (for φ^{*,j}). + * @param t The source time index (φᵗ). + * @param vec The vector v to scale. + * @return The scaled vector γₘ · v, or zero if m is out of bounds. + */ + template + VectorType + apply_derivative( + int const &j, + int const &t, + VectorType const &vec + ) const { + + int const r = std::min(j, r_order); // Effective order. + int const m = (j - 1) - t; // Index difference m = j−1−t + + VectorType result(vec); // Copy the input vector v. + + if (m < 0 || m >= r) { + result = 0.0; // Zero contribution if out of bounds. + } else { + auto gammas = this->compute_gammas(r); + result *= gammas[m]; // Scale by γₘ. + } + + return result; + + } + +private: + + /// @brief The finite element space in which the extrapolated fields are defined. + mfem::ParFiniteElementSpace *space_; + + /** + * @brief Compute extrapolation weights γₘ for φ^{*,k+1} given an extrapolation order. + * + * These weights define the extrapolated value: + * φ^{*,k+1} = ∑ₘ γₘ · φ^{k−m}, for m = 0..order−1. + * + * @param order The extrapolation order (must be 0, 1, or 2). + * @return A vector of γₘ weights indexed by m. + */ + std::vector + compute_gammas( + int order + ) const { + + SLIC_ERROR_ROOT_IF( + order < 0 || order > 2, + axom::fmt::format("Unsupported extrapolation order = {}. Only orders 0, 1, and 2 are supported.", order) + ); + + switch (order) { + case 0: + return {}; // Zero extrapolation → φ^{*,k+1} = 0 + case 1: + return {1.0}; // Constant extrapolation → φ^{*,k+1} = φ^k + case 2: + return {2.0, -1.0}; // Linear extrapolation → φ^{*,k+1} = 2φ^k − φ^{k−1} + default: + // Should be unreachable due to earlier guard. + SLIC_ERROR_ROOT("Unreachable case in switch for extrapolation weights."); + return {}; + } + + } + +}; // struct ExtrapolationOperator + +} // namespace serac diff --git a/src/serac/numerics/fixed_step_bdf_operator.hpp b/src/serac/numerics/fixed_step_bdf_operator.hpp new file mode 100644 index 000000000..5215c3511 --- /dev/null +++ b/src/serac/numerics/fixed_step_bdf_operator.hpp @@ -0,0 +1,216 @@ +// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and +// other Serac Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/** + * @file fixed_step_bdf_operator.hpp + * + * @brief Implements a Backward Differentiation Formula (BDF) time derivative + * operator with fixed timestep and dynamic order selection. + */ + +#pragma once + +#include "serac/infrastructure/logger.hpp" + +namespace serac { + +/** + * @brief Implements a fixed time step Backward Differentiation Formula (BDF) operator. + * + * This class supports BDF1 through BDF6 schemes with constant timestep. It computes + * the BDF coefficients dynamically based on how many previous states are available. + * + * @tparam s_order The maximum order of the BDF method (e.g., 1 to 6). + */ +template +struct FixedStepBDFOperator { + + /** + * @brief Computes the BDF time derivative: Dφ^{k+1} = (β₀ * φ^{k+1} + weighted_sum_previous) / Δt + * + * The BDF method approximates the time derivative at time step k+1 using a weighted + * combination of current and past states. The `weighted_sum_previous` argument is assumed + * to already include the correct signs for β₁ through β_q. + * + * As a result, this function uses a '+' to combine β₀ * φ^{k+1} with the weighted sum. + * + * @tparam T1 Type of the current state φ^{k+1}, includes AD support (e.g., dual numbers). + * @tparam T2 Type of the weighted sum of previous states (typically plain). + * + * @param phi_k_plus_1 The current state at time step k+1. + * @param weighted_sum_previous The precomputed sum β₁ * φ^k + β₂ * φ^{k−1} + … (signs included). + * @param dt The time step size Δt. + * @param k_plus_1 Current time step index (i.e., solving for time step \( k + 1 \)) + * + * @return The time derivative Dφ^{k+1} + */ + template + auto + time_derivative( + T1 phi_k_plus_1, // e.g., serac::tensor, 2> + T2 weighted_sum_previous, // e.g., serac::tensor + double dt, + int k_plus_1 + ) const { + + int available_order = std::min(k_plus_1, s_order); + auto beta = compute_beta(available_order); + + SLIC_ERROR_ROOT_IF( + static_cast(beta.size()) < 1, + "Insufficient BDF coefficients computed. Cannot evaluate time derivative." + ); + + // Compute the full BDF expression: (β₀ * φ^{k+1} + Σ_j β_j * φ^{k-j}) / Δt + return ((beta[0] * phi_k_plus_1) + weighted_sum_previous) / dt; + + } + + /** + * @brief Computes the weighted sum of previous states using dynamic BDF coefficients. + * + * @param previous_states A container of previous states order from newest (`u^k`) to oldest (`u^{k-q+1}`). + * @param k_plus_1 Current time step index (i.e., solving for time step \( k + 1 \)) + * + * @return The weighted sum of previous values. + */ + ::serac::FiniteElementState + weighted_sum( + std::vector<::serac::FiniteElementState> const &previous_states, + int k_plus_1 + ) const { + + SLIC_ERROR_ROOT_IF( + k_plus_1 == 0 || previous_states.empty(), + "No previous states available to compute weighted sum." + ); + + int available_order = std::min(k_plus_1, s_order); + + SLIC_ERROR_ROOT_IF( + static_cast(previous_states.size()) != available_order, + ::axom::fmt::format("Expected {} previous states, but got {}.", available_order, previous_states.size()) + ); + + auto const beta = compute_beta(available_order); + + ::serac::FiniteElementState result = previous_states[0]; + result *= beta[1]; // Apply first beta coefficient (skip beta[0], which is for u_{k+1}) + + for (int i = 1; i < available_order; ++i) { + result.Add(beta[i + 1], previous_states[i]); + } + + return result; + + } + + /** + * @brief Retrieves the BDF coefficient β_j for a given cycle and index. + * + * @param k_plus_1 Current time step index (i.e., solving for time step \( k + 1 \)) + * @param j Index of the coefficient (0 = next-step, 1 = previous, etc.). + * @return The β_j coefficient for the current order. + */ + double + compute_beta( + int k_plus_1, + int j + ) const { + + SLIC_ERROR_ROOT_IF(k_plus_1 <= 0, "Cannot compute BDF coefficient: cycle must be greater than zero."); + + int available_order = std::min(k_plus_1, s_order); + auto betas = compute_beta(available_order); + + SLIC_ERROR_ROOT_IF( + j < 0 || j >= static_cast(betas.size()), + ::axom::fmt::format("Index {} out of range for BDF coefficients (available: {}).", j, betas.size()) + ); + + return betas[j]; + + } + + /** + * @brief Applies the derivative ∂(Dφ^j)/∂φ^t to a vector. + * + * This computes the action of the BDF derivative’s Jacobian with respect to φ^t: + * ∂/∂φ^t [ (β₀ * φ^j + ∑ₘ βₘ * φ^{j−m}) / Δt ] = (βₘ / Δt) · I + * where m = j − t and t ∈ {j, j−1, ..., j−(r−1)}. + * + * If t is outside this range, the result is zero. + * + * @tparam VectorType A vector-like type supporting assignment and scalar multiplication. + * @param j The current step index (computing ∂/∂φ^t of Dφ^j). + * @param t The time index of the derivative target φ^t. + * @param vec The vector v to scale. + * @param dt The time step Δt. + * @return A scaled copy of vec, or zero if out of valid range. + */ + template + VectorType + apply_derivative( + int j, + int t, + VectorType const &vec + ) const { + + int r = std::min(j, s_order); // Effective BDF order. + int m = j - t; + + VectorType result(vec); // Copy input vector. + + if (m < 0 || m >= r + 1) { + result = 0.0; + } else { + auto beta = compute_beta(r); // Get β₀ through β_r. + result *= beta[m]; + } + + return result; + + } + + /** + * @brief Returns the BDF order actually in use. + */ + static int constexpr order() { return s_order; } + +private: + + /** + * @brief Internal function to compute all β_j coefficients for a given order. + * + * @param order The order of the BDF method (≤ s_order). + * @return Vector of BDF coefficients {β₀, β₁, ..., β_order}. + */ + std::vector + compute_beta( + int order + ) const { + + SLIC_ERROR_ROOT_IF( + order < 1 || order > 6, + ::axom::fmt::format("Unsupported BDF order: {} (valid range is 1 to 6).", order) + ); + + switch (order) { + case 1: return {1.0, -1.0}; + case 2: return {3.0 / 2.0, -2.0, 0.5}; + case 3: return {11.0 / 6.0, -3.0, 1.5, -1.0 / 3.0}; + case 4: return {25.0 / 12.0, -4.0, 3.0, -4.0 / 3.0, 0.25}; + case 5: return {137.0 / 60.0, -5.0, 5.0, -10.0 / 3.0, 1.25, -0.2}; + case 6: return {147.0 / 60.0, -6.0, 7.5, -20.0 / 3.0, 3.75, -1.2, 1.0 / 6.0}; + default: + return {}; + } + + } + +}; // struct FixedStepBDFOperator + +} // namespace serac diff --git a/src/serac/numerics/functional/tensor.hpp b/src/serac/numerics/functional/tensor.hpp index 83e5a6b41..075720dc2 100644 --- a/src/serac/numerics/functional/tensor.hpp +++ b/src/serac/numerics/functional/tensor.hpp @@ -85,6 +85,14 @@ struct tensor { return data[0]; } + SERAC_HOST_DEVICE constexpr tensor& operator=(const T& scalar) + { + for (int i = 0; i < m; ++i) { + data[i] = scalar; + } + return *this; + } + T data[m]; }; /// @endcond @@ -1338,7 +1346,7 @@ SERAC_HOST_DEVICE constexpr auto detApIm1(const tensor& A) // clang-format off // equivalent to tr(A) + I2(A) + det(A) - return A(0, 0) + A(1, 1) + A(2, 2) + return A(0, 0) + A(1, 1) + A(2, 2) - A(0, 1) * A(1, 0) * (1 + A(2, 2)) + A(0, 0) * A(1, 1) * (1 + A(2, 2)) - A(0, 2) * A(2, 0) * (1 + A(1, 1)) @@ -1946,6 +1954,38 @@ SERAC_HOST_DEVICE constexpr int size(const double&) { return 1; } /// @overload SERAC_HOST_DEVICE constexpr int size(zero) { return 0; } +/** + * @brief Returns the rank (number of dimensions) of a tensor. + * + * @tparam T The datatype stored in the tensor. + * @tparam n The extents of each dimension. + * + * @return The number of dimensions (rank) of the tensor. + */ +template +SERAC_HOST_DEVICE int constexpr +rank( + ::serac::tensor const & +) { + return sizeof...(n); +} + +/** + * @brief Returns the shape (extents of each dimension) of a tensor. + * + * @tparam T The datatype stored in the tensor. + * @tparam n The extents of each dimension. + * + * @return A ::std::array of extents representing the tensor shape. + */ +template +SERAC_HOST_DEVICE ::std::array constexpr +shape( + ::serac::tensor const & +) { + return {n...}; +} + /** * @brief a function for querying the ith dimension of a tensor * @@ -2089,7 +2129,7 @@ inline mat < 2, 2 > look_at(const vec < 2 > & direction) { inline mat < 3, 3 > R3_basis(const vec3 & n) { float sign = (n[2] >= 0.0f) ? 1.0f : -1.0f; - float a = -1.0f / (sign + n[2]); + float a = -1.0f / (sign + n[2]); float b = n[0] * n[1] * a; return mat < 3, 3 >{ diff --git a/src/serac/physics/base_physics.cpp b/src/serac/physics/base_physics.cpp index 5af5d0145..982d91712 100644 --- a/src/serac/physics/base_physics.cpp +++ b/src/serac/physics/base_physics.cpp @@ -386,6 +386,88 @@ std::unordered_map BasePhysics::getCheckpointed return previous_states_map; } +std::vector +BasePhysics::getPreviousStates( + int const &cycle, + int const &num_requested_states, + std::string const &field_name +) { + + int const num_available = std::min(num_requested_states, cycle); + std::vector<::serac::FiniteElementState> previous_states; + + for (int i = 0; i < num_available; ++i) { + + int checkpoint_step = cycle - 1 - i; + + SLIC_ERROR_ROOT_IF( + checkpoint_step < 0, + axom::fmt::format("Invalid checkpoint step: {} (cycle = {})", checkpoint_step, cycle) + ); + + // First check: Does the field exist in the overall checkpoint map? + SLIC_ERROR_ROOT_IF( + checkpoint_states_.find(field_name) == checkpoint_states_.end(), + axom::fmt::format("Field '{}' was never registered in the checkpoint history.", field_name) + ); + + // Second check: Is the checkpoint step valid for this field? + SLIC_ERROR_ROOT_IF( + checkpoint_step >= static_cast(checkpoint_states_.at(field_name).size()), + axom::fmt::format("Checkpoint step {} exceeds stored checkpoints for field '{}'.", checkpoint_step, field_name) + ); + + // Now we can safely fetch that particular timestep and check the field within it. + auto const &checkpointed = this->getCheckpointedStates(checkpoint_step); + + previous_states.push_back(checkpointed.at(field_name)); + + } + + return previous_states; + +} + +double +BasePhysics::getAbsoluteChangeState( + int const &cycle, + std::string const &state_name +) { + + SLIC_ERROR_ROOT_IF( + cycle < 1, + axom::fmt::format( + "Cannot compute absolute change for state '{}': invalid cycle {} (must be >= 1).", + state_name, + cycle + ) + ); + + auto current_state = this->getCheckpointedStates(cycle).at(state_name); + auto previous_state = this->getCheckpointedStates(cycle - 1).at(state_name); + + auto diff(current_state); + diff.Add(-1.0, previous_state); + + return norm(diff); + +} + +double +BasePhysics::getRelativeChangeState( + int const &cycle, + std::string const &state_name +) { + + double norm_diff = this->getAbsoluteChangeState(cycle, state_name); + + auto current_state = this->getCheckpointedStates(cycle).at(state_name); + double norm_curr = norm(current_state); + + return norm_diff / (norm_curr + 1.0e-12); // Prevent division by zero. + +} + double BasePhysics::getCheckpointedTimestep(int cycle) const { SLIC_ERROR_ROOT_IF(cycle < 0, axom::fmt::format("Negative cycle number requested for physics module {}.", name_)); @@ -395,6 +477,20 @@ double BasePhysics::getCheckpointedTimestep(int cycle) const return cycle < static_cast(timesteps_.size()) ? timesteps_[static_cast(cycle)] : 0.0; } +double BasePhysics::getCheckpointedTime(int cycle) const +{ + SLIC_ERROR_ROOT_IF(cycle < 0, axom::fmt::format("Negative cycle number requested for physics module {}.", name_)); + SLIC_ERROR_ROOT_IF(cycle > static_cast(timesteps_.size()), + axom::fmt::format("Time for cycle {} requested, but physics module has only reached cycle {}.", + cycle, timesteps_.size())); + + double time = 0.0; + for (int i = 0; i < cycle; ++i) { + time += timesteps_[static_cast(i)]; + } + return time; +} + namespace detail { std::string addPrefix(const std::string& prefix, const std::string& target) { diff --git a/src/serac/physics/base_physics.hpp b/src/serac/physics/base_physics.hpp index 138e0fdcc..a9105051e 100644 --- a/src/serac/physics/base_physics.hpp +++ b/src/serac/physics/base_physics.hpp @@ -446,10 +446,51 @@ class BasePhysics { /** * @brief Get a timestep increment which has been previously checkpointed at the give cycle * @param cycle The previous 'timestep' number where the timestep increment is requested - * @return The timestep increment + * @return The timestep increment. */ virtual double getCheckpointedTimestep(int cycle) const; + /** + * @brief Get the cumulative time up to the given cycle by summing previously checkpointed timestep increments. + * @param cycle The previous 'timestep' number where the simulation time is requested + * @return The total simulation time at the start of the given cycle. + */ + virtual double getCheckpointedTime(int cycle) const; + + /** + * @brief Computes the absolute change (L2 norm of the difference) between two successive states. + * + * This function compares checkpointed states at cycles `k` and `k-1` for a given field. + * + * @param cycle The solution cycle `k` (must be >= 1). + * @param state_name The name of the checkpointed state field to compare. + * + * @return The absolute L2 norm of the difference: || u^k - u^{k-1} ||. + */ + double + getAbsoluteChangeState( + int const &cycle, + std::string const &state_name + ); + + /** + * @brief Computes the relative change in a given state variable between successive cycles. + * + * This returns the L2 norm of the difference divided by the L2 norm of the current state: + * + * || u^k - u^{k-1} || / (|| u^k || + ε) + * + * @param cycle The solution cycle `k` (must be >= 1). + * @param state_name The name of the checkpointed state field to compare. + * + * @return The relative L2 change between `state_name` at cycles `k` and `k-1`. + */ + double + getRelativeChangeState( + int const &cycle, + std::string const &state_name + ); + /** * @brief Initializes the Sidre structure for simulation summary data * @@ -530,6 +571,27 @@ class BasePhysics { */ std::unordered_map getCheckpointedStates(int cycle); + /** + * @brief Retrieves a set of previous state values for time integration schemes + * (e.g., BDF, Extrapolation). + * + * This version pulls data directly from the solver's checkpointing system via + * `getCheckpointedStates(cycle)`. The returned states are ordered from newest + * (`u^k`) to oldest (`u^{k-q+1}`). + * + * @param cycle The current timestep index `k+1` when solving for `u^{k+1}`. + * @param num_requested_states Number of past states required (e.g., 2 for BDF2, 1 for linear extrapolation). + * @param field_name The key identifying the field in checkpointed states (e.g., "velocity", "pressure"). + * + * @return A vector of previous states, ordered from latest to oldest. + */ + std::vector + getPreviousStates( + int const &cycle, + int const &num_requested_states, + std::string const &field_name + ); + /// @brief Name of the physics module std::string name_ = {}; diff --git a/src/serac/physics/materials/thermal_material.hpp b/src/serac/physics/materials/thermal_material.hpp index 98606c7ba..3cf1f864d 100644 --- a/src/serac/physics/materials/thermal_material.hpp +++ b/src/serac/physics/materials/thermal_material.hpp @@ -65,6 +65,67 @@ struct LinearIsotropicConductor { double conductivity_; }; +struct ConstantIsotropicDiffuser { + + /** + * @brief Construct a new Linear Isotropic Diffuser object + * + * @param density Density of the carrier fluid/material (mass/volume). + * @param diffusivity Mass diffusivity (length^2 / time). + */ + ConstantIsotropicDiffuser( + double density = 1.0, + double diffusivity = 1.0 + ) : density_(density) + , diffusivity_(diffusivity) { + + SLIC_ERROR_ROOT_IF( + density_ < 0.0, + "Density must be positive in the linear isotropic diffuser model." + ); + + SLIC_ERROR_ROOT_IF( + diffusivity_ < 0.0, + "Diffusivity must be positive in the linear isotropic diffuser model." + ); + + } + + /** + * @brief Material response call for a linear isotropic diffusion material. + * + * @tparam T1 Spatial position type. + * @tparam T2 Scalar field type (e.g., species mass fraction Y). + * @tparam T3 Scalar gradient type. + + * @param[in] gradient_y Gradient of the transported scalar Y. + + * @return The material response (tuple of capacity and diffusive flux). + */ + template + SERAC_HOST_DEVICE auto + operator()( + T1 const& /* x */, + T2 const& /* Y */, + T3 const& gradient_Y + ) const { + + // Capacity = ρ * C (with C = 1 for species) + // Flux = -ρ D ∇Y + return serac::tuple{density_, 1.0, -1.0 * density_ * diffusivity_ * gradient_Y}; + + } + + private: + + /// Constant density. + double density_; + + /// Constant isotropic mass diffusivity. + double diffusivity_; + +}; + /// Nonlinear isotropic heat transfer material model struct IsotropicConductorWithLinearConductivityVsTemperature { /** diff --git a/src/serac/physics/state/finite_element_state.hpp b/src/serac/physics/state/finite_element_state.hpp index 0132ee0e2..0780091f4 100644 --- a/src/serac/physics/state/finite_element_state.hpp +++ b/src/serac/physics/state/finite_element_state.hpp @@ -82,7 +82,10 @@ using first_argument = std::decay_t auto evaluateTensorFunctionOnMfemVector(const mfem::Vector& X_mfem, Callable&& f) { - first_argument X; + using Arg = first_argument; + using Bare = std::decay_t; // strip const& + + Bare X{}; SLIC_ERROR_IF(X_mfem.Size() != size(X), "Size of tensor in callable does not match spatial dimension of MFEM Vector."); for (int i = 0; i < X_mfem.Size(); i++) X[i] = X_mfem[i]; diff --git a/src/serac/physics/state/finite_element_vector.cpp b/src/serac/physics/state/finite_element_vector.cpp index 917e37559..f0dbed151 100644 --- a/src/serac/physics/state/finite_element_vector.cpp +++ b/src/serac/physics/state/finite_element_vector.cpp @@ -149,4 +149,17 @@ double innerProduct(const FiniteElementVector& v1, const FiniteElementVector& v2 return global_ip; } +double norml2(const mfem::HypreParVector& fe_vector) +{ + // GlobalVector() allocates a new mfem::Vector on the heap containing all entries + // from the parallel vector, and transfers ownership of it to the caller. To + // avoid confusion and ensure this memory is released automatically when we exit + // this scope, we wrap it in a std::unique_ptr. + std::unique_ptr global_vec(fe_vector.GlobalVector()); + + // Compute and return the Euclidean L2 norm of the fully assembled global vector. + // Norml2() computes sqrt(sum_i v_i^2) over all entries. + return global_vec->Norml2(); +} + } // namespace serac diff --git a/src/serac/physics/state/finite_element_vector.hpp b/src/serac/physics/state/finite_element_vector.hpp index 6421b6d35..cd164f725 100644 --- a/src/serac/physics/state/finite_element_vector.hpp +++ b/src/serac/physics/state/finite_element_vector.hpp @@ -235,4 +235,16 @@ double min(const FiniteElementVector& fe_vector); */ double innerProduct(const FiniteElementVector& vec1, const FiniteElementVector& vec2); +/** + * @brief Compute the global L2 norm of a distributed finite element vector. + * + * This function collects the distributed entries of an mfem::HypreParVector + * into a single global mfem::Vector on the root processor, then computes its + * Euclidean (L2) norm. + * + * @param fe_vector The distributed finite element vector (HypreParVector). + * @return The global L2 norm (sqrt of sum of squares) of the vector. + */ +double norml2(const mfem::HypreParVector& fe_vector); + } // namespace serac diff --git a/src/serac/physics/state/state_manager.cpp b/src/serac/physics/state/state_manager.cpp index e3a400d7d..85a5eeeb5 100644 --- a/src/serac/physics/state/state_manager.cpp +++ b/src/serac/physics/state/state_manager.cpp @@ -195,7 +195,7 @@ void StateManager::save(const double t, const int cycle, const std::string& mesh auto& datacoll = datacolls_.at(mesh_tag); std::string file_path = axom::utilities::filesystem::joinPath(datacoll.GetPrefixPath(), datacoll.GetCollectionName()); SLIC_INFO_ROOT( - axom::fmt::format("Saving data collection at time: '{}' and cycle: '{}' to path: '{}'", t, cycle, file_path)); + axom::fmt::format("Saving data collection at time: '{0:1.12e}' and cycle: '{1}' to path: '{2}'", t, cycle, file_path)); datacoll.SetTime(t); datacoll.SetCycle(cycle); diff --git a/src/serac/physics/state/state_manager.hpp b/src/serac/physics/state/state_manager.hpp index 2ce71fc38..2913bbeed 100644 --- a/src/serac/physics/state/state_manager.hpp +++ b/src/serac/physics/state/state_manager.hpp @@ -275,6 +275,42 @@ class StateManager { */ static FiniteElementDual newDual(const mfem::ParFiniteElementSpace& space, const std::string& dual_name); + /** + * @brief Factory method for creating an array of new FEDual objects. + * + * @tparam N The number of duals to create. + * @tparam FunctionSpace The function space (e.g. H1<1>) to build the finite element duals on. + * + * @param space The function space used to build the finite element duals. + * @param full_prefix A common string prefix for naming the duals (each dual will be named as + * "_", where runs from 0 to N-1). + * @param mesh_tag The tag for the stored mesh used to construct the finite element duals. + * + * @return An std::array of N constructed finite element duals. + * + * @see StateManager::newDual + * @note If this is a restart then the options (except for the names) will be ignored. + */ + template + static std::array + newDualArray( + FunctionSpace space, + std::string const &full_prefix, + std::string const &mesh_tag + ) { + + return [&](std::index_sequence) { + return std::array{ + StateManager::newDual( + space, + full_prefix + "_" + std::to_string(I), // Just append index. + mesh_tag + )... + }; + }(std::make_index_sequence{}); + + } + /** * @brief Store a pre-constructed finite element dual in the state manager *