feat(equations): added mass normalization

This equation weakly constrains the total mass of the system to some prescribed mass by integrating the density over the volume

src/include/serif/equations/mass_normalization.hpp

@@ -0,0 +1,110 @@
+#pragma once
+
+/**
+ * @file mass_normalization.hpp
+ * @author Emily M. Boudreaux
+ *
+ * @brief File containing the mass normalization equation, which enforces a global constraint on the total mass of a system.
+ */
+
+#include <cmath>
+#include <string_view>
+
+#include "serif/equations/contributions.hpp"
+#include "serif/equations/point.hpp"
+#include "serif/utils/error/errors.hpp"
+
+namespace serif::equations {
+
+    /**
+     * @brief Mass normalization equation, which enforces a global constraint on the total mass of a system.
+     *
+     *
+     * @note In the barotropic case the bernoulli constant as the control variable here. This is a "clever" choice in this case; however, more broadly most
+     * of this equation is general (its only that in we choose to inject the residual
+     * into the mass_normalization member of the barotropic_constant block). That could very easily become a more
+     * general control variable without changing the rest of this
+     *
+     *
+     * Mass normalization, being an integral constraint must make use of the ::residual_offset() method to properly find the residual. See AngularMomentumConstraint
+     * for a more detailed discussion of this. Suffice it to say that these are point-wise integrands so if we were to subtract the
+     * target mass from every point we would wildly over count the target mass. Rather, we simply return the density everywhere as the
+     * volume integrated density will be the mass. Then the ::residual_offset() method returns the negative of the target mass
+     * (the backend automatically adds anything from this function, if it exists, to the integrated residual to get the final residual).
+     *
+     * @contribution Residuals to the mass_normalization term in the barotropic_constant block. This is effectively a global constraint on the total mass of the system.
+     */
+    class MassNormalization final : public Contribution<discretization::forms::Operand<discretization::blocks::barotropic_constant::mass_normalization>> {
+    public:
+        using Density = discretization::blocks::density::mass;          ///< Alias for mass density
+        using ResidualQuantity = dimensions::Mass;                      ///< Alias for the dimension of the residual
+        using DensityValue = discretization::forms::Operand<Density>;   ///< Alias for the operand type of the mass density
+        using Operands = utils::types::TypeList<DensityValue>;          ///< The list of operands for the mass normalization equation.
+
+        static constexpr std::string_view name = "mass normalization";  ///< Compile time name for the mass normalization equation, used for debugging and logging.
+
+        /**
+         * @brief Construct a mass normalization equation.
+         *
+         * @param target_mass The target mass for the system, which must be finite and positive.
+         *
+         * @throws std::invalid_argument If the target mass is not finite or not positive.
+         */
+        explicit MassNormalization(const double target_mass)
+        : m_target_mass(target_mass) {
+            using namespace utils::errors;
+
+            serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(std::isfinite(target_mass) && target_mass > 0.0, "Target mass must be finite and positive.");
+        }
+
+        /**
+         * @brief Evaluate the mass normalization equation point-wise integrand and its partial derivatives.
+         *
+         * The point-wise integrand is calculated as \f$ \rho \f$, where \f$ \rho \f$ is the model density. The partial derivative is calculated as:
+         * - \f$ \frac{\partial r}{\partial \rho} = 1 \f$
+         *
+         * @note As the class description notes, the residual is not calculated as \f$ \rho - M_{\text{target}} \f$ because this would over count the target mass. Instead we simply return the density and then use the ::residual_offset() method to subtract the target mass from the integrated residual.
+         *
+         * @param point The point at which to evaluate the equation.
+         * @returns The point-wise linearization of the mass normalization equation, which contains the residual and the partial derivative with respect to the model density.
+         */
+        [[nodiscard]] static PointLinearization<Operands> evaluate(const EquationPoint<Operands>& point) {
+            PointLinearization<Operands> result;
+            result.residual = point.fields.get<Density>();
+            result.field_derivatives.get<Density>() = 1.0;
+            return result;
+        }
+
+        /**
+         * @brief Return the additional quadrature order required for the operands of this equation.
+         *
+         * @param operand The operand type for which to determine the additional quadrature order. Note that this is only used to select the correct overload of this function. This particular version only has  one overload; therefore this is really just their to maintain symmetry with other
+         * equation implementations.
+         * @param position_order The order of the position basis functions.
+         * @return The additional quadrature order required for the given operand type (**always 0 in this case**).
+         */
+        [[nodiscard]] static int additional_order(DensityValue /*order operand type*/, int) {
+            return 0;
+        }
+
+        /**
+         * @brief Return the residual offset for the mass normalization equation.
+         *
+         * The backend integrates the point-wise residual integrand over the domain. The global constraint
+         * is satisfied when the total mass equals the target value. Therefore we subtract the target
+         * from the integrated residual to get the actual residual.
+         *
+         * @note The backend automatically looks for this method and calls it if it exists. <br><br>
+         *
+         * @note The backend actually adds whatever is returned from this function to the integrated residual to get the final residual. Therefore we return the negative of the target mass here.
+         *
+         * @return The negative of the target mass, $-M_{\text{target}}$.
+         */
+        [[nodiscard]] double residual_offset() const {
+            return -m_target_mass;
+        }
+
+    private:
+        double m_target_mass;  ///< The target mass for the system, which must be finite and positive.
+    };
+} // namespace serif::equations