manning_roughness_1d.h

Go to the documentation of this file.

//
// Created by julian on 1/29/25.
//
 
#ifndef D3F_MANNING_ROUGHNESS_1D_H
#define D3F_MANNING_ROUGHNESS_1D_H
#include "lib_disc/spatial_disc/user_data/linker/linker.h"
 
namespace ug {namespace d3f{
template<int dim>
class ManningRoughness1d : public StdDataLinker<ManningRoughness1d<dim>, number, dim> {
 
    // --- Type Definitions ---
    typedef CplUserData<number, dim>             TNumberData;
    typedef CplUserData<MathVector<dim>, dim>    TVectorData;
    typedef DependentUserData<number, dim>       TDNumberData;
    typedef DependentUserData<MathVector<dim>, dim> TDVectorData;
 
public:
    // Identifiers for the input fields
    enum inputs { _h_ = 0, _v_, _n_ };
 
    ManningRoughness1d() {
        this->set_num_input(3);
    }
 
    // --- Setters ---
    void set_hmin(number h)      { m_hmin = h; } // CRITICAL FIX: Removed stray 'return;'
    void set_gravity(number gin) { m_g = gin; }
 
    void set_manning_coefficient(SmartPtr<TNumberData> n) {
        m_n = n;
        m_dn = n.template cast_dynamic<TDNumberData>();
        this->set_input(_n_, n, n);
    }
 
    void set_h(SmartPtr<TNumberData> h) {
        m_h = h;
        m_dh = h.template cast_dynamic<TDNumberData>();
        this->set_input(_h_, h, h);
    }
 
    void set_v(SmartPtr<TNumberData> v) {
        m_v = v;
        m_dv = v.template cast_dynamic<TDNumberData>();
        this->set_input(_v_, v, v);
    }
 
    // --- Evaluation Methods ---
 
    inline void evaluate(number& value, const MathVector<dim>& globIP, number time, int si) const {
        UG_THROW("ManningRoughness1d: Not implemented for single global positions.\n");
    }
 
    template <int refDim>
    inline void evaluate(number vValue[],
                         const MathVector<dim> vGlobIP[],
                         number time, int si,
                         GridObject* elem,
                         const MathVector<dim> vCornerCoords[],
                         const MathVector<refDim> vLocIP[],
                         const size_t nip,
                         LocalVector* u,
                         const MathMatrix<refDim, dim>* vJT = NULL) const
    {
        // Note: Using VLAs (Variable Length Arrays) here. standard in UG4, but be aware it's a C99/compiler extension.
        number vVel[nip];
        number vh[nip];
        number vn[nip];
 
        // Evaluate inputs at the integration points
        (*m_h)(&vh[0], vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
        (*m_v)(&vVel[0], vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
        (*m_n)(&vn[0], vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
 
        for (size_t ip = 0; ip < nip; ++ip) {
            // Cutoff to prevent division by zero or oscillating signs at near-zero velocities
            if (std::fabs(vVel[ip]) < m_vmin) {
                vValue[ip] = 0.0;
                continue;
            }
 
            // Calculate Reaction Rate R = (g * n^2 * |v|) / h^(4/3)
            vValue[ip] = m_g * vn[ip] * vn[ip] * std::fabs(vVel[ip]);
            vValue[ip] /= std::pow(std::max(vh[ip], m_hmin), (4.0 / 3.0));
        }
    }
 
    template <int refDim>
    void eval_and_deriv(number vValue[],
                        const MathVector<dim> vGlobIP[],
                        number time, int si,
                        GridObject* elem,
                        const MathVector<dim> vCornerCoords[],
                        const MathVector<refDim> vLocIP[],
                        const size_t nip,
                        LocalVector* u,
                        bool bDeriv,
                        int s,
                        std::vector<std::vector<number>> vvvDeriv[],
                        const MathMatrix<refDim, dim>* vJT = NULL) const
    {
        // 1. Evaluate the base value
        evaluate<refDim>(vValue, vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
 
        if (!bDeriv) return;
        if (this->zero_derivative()) return; // CRITICAL FIX: Hoisted out of the integration point loop
 
        number vVel[nip];
        number vh[nip];
        number vn[nip];
 
        // Re-evaluate inputs to use for derivative calculations
        (*m_h)(&vh[0], vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
        (*m_v)(&vVel[0], vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
        (*m_n)(&vn[0], vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
 
        this->set_zero(vvvDeriv, nip);
 
        for (size_t ip = 0; ip < nip; ++ip) {
            // Calculate partial derivatives analytically
            number dh = evaluate_dh(vVel[ip], vn[ip], vh[ip], m_g);
            number dv = evaluate_dv(vVel[ip], vn[ip], vh[ip], m_g);
            number dn = evaluate_dn(vVel[ip], vn[ip], vh[ip], m_g);
 
            // --- Derivative w.r.t Depth (h) ---
            size_t sid = this->series_id(_h_, s);
            if (m_dh.valid() && !m_dh->zero_derivative()) {
                for (size_t fct = 0; fct < m_dh->num_fct(); fct++) {
                    const number* vDHeight = m_dh->deriv(sid, ip, fct);
                    const size_t commonFct = this->input_common_fct(_h_, fct);
 
                    if (this->num_sh(commonFct) == nip) {
                        vvvDeriv[ip][commonFct][ip] += dh * vDHeight[ip]; // FV1 mass lumping
                    } else {
                        for (size_t sh = 0; sh < this->num_sh(commonFct); ++sh) {
                            vvvDeriv[ip][commonFct][sh] += dh * vDHeight[sh];
                        }
                    }
                }
            }
 
            // --- Derivative w.r.t Velocity (v) ---
            sid = this->series_id(_v_, s);
            if (m_dv.valid() && !m_dv->zero_derivative()) {
                for (size_t fct = 0; fct < m_dv->num_fct(); fct++) {
                    const number* vDVel = m_dv->deriv(sid, ip, fct);
                    const size_t commonFct = this->input_common_fct(_v_, fct);
 
                    if (this->num_sh(commonFct) == nip) {
                        vvvDeriv[ip][commonFct][ip] += dv * vDVel[ip];
                    } else {
                        for (size_t sh = 0; sh < this->num_sh(commonFct); ++sh) {
                            vvvDeriv[ip][commonFct][sh] += dv * vDVel[sh];
                        }
                    }
                }
            }
 
            // --- Derivative w.r.t Manning's n ---
            sid = this->series_id(_n_, s);
            if (m_n.valid() && !m_n->zero_derivative()) {
                for (size_t fct = 0; fct < m_dn->num_fct(); fct++) {
                    const number* vDn = m_dn->deriv(sid, ip, fct);
                    const size_t commonFct = this->input_common_fct(_n_, fct);
 
                    if (this->num_sh(commonFct) == nip) {
                        vvvDeriv[ip][commonFct][ip] += dn * vDn[ip];
                    } else {
                        for (size_t sh = 0; sh < this->num_sh(commonFct); ++sh) {
                            vvvDeriv[ip][commonFct][sh] += dn * vDn[sh];
                        }
                    }
                }
            }
        }
    }
 
    bool requires_grid_fct() const { return true; }
 
    // --- Analytical Derivatives ---
 
    // Partial derivative w.r.t velocity: dR/dv
    number evaluate_dv(const number v, const number n, const number h, const number g) const {
        // v / fabs(v + SMALL_SQ) acts as a smoothed sign(v) function
        return g * n * n * v / (std::fabs(v + SMALL_SQ) * std::pow(std::max(h, m_hmin), (4.0 / 3.0)));
    }
 
    // Partial derivative w.r.t depth: dR/dh
    number evaluate_dh(const number v, const number n, const number h, const number g) const {
        return (-4.0 / 3.0) * g * n * n * std::fabs(v) / std::pow(std::max(h, m_hmin), (7.0 / 3.0));
    }
 
    // Partial derivative w.r.t Manning's n: dR/dn
    number evaluate_dn(const number v, const number n, const number h, const number g) const {
        return (2.0 * g * n * std::fabs(v)) / std::pow(std::max(h, m_hmin), (4.0 / 3.0));
    }
 
private:
    SmartPtr<TNumberData>  m_v;
    SmartPtr<TDNumberData> m_dv;
    SmartPtr<TNumberData>  m_h;
    SmartPtr<TDNumberData> m_dh;
    SmartPtr<TNumberData>  m_n;
    SmartPtr<TDNumberData> m_dn;
 
    // Member variables standardized with 'm_' prefix
    number m_vmin = 1e-6;
    number m_hmin = 1e-6;
    number m_g    = 9.81 * 60.0 * 60.0; // Gravity, scaled likely for time units (e.g., hours)
};
 
 
}
}// namespace ug
#endif //MANNING_ROUGHNESS_1D_H