st_venant_elem_disc.h

Go to the documentation of this file.

//
// Created by julian on 5/15/25.
//
 
#ifndef D3F_ST_VENANT_ELEM_DISC_H
#define D3F_ST_VENANT_ELEM_DISC_H
#include "manning_roughness_1d.h"
#include "sohle.h"
#include "../../../ConvectionDiffusion/fv1/convection_diffusion_fv1.h"
#include "lib_disc/spatial_disc/user_data/linker/inverse_linker.h"
#include "lib_disc/spatial_disc/user_data/linker/scale_add_linker.h"
#include "lib_disc/spatial_disc/disc_util/conv_shape.h"
#include "lib_disc/spatial_disc/user_data/common_user_data/element_orientation_data.h"
#include "st_venant_wall_bnd.h"
 
namespace ug {namespace d3f{
 
template <typename TDomain>
class StVenantElemDisc : public ug::ConvectionDiffusionPlugin::ConvectionDiffusionFV1<TDomain> {
public:
    static const int dim = TDomain::dim;
 
    // --- Type Definitions ---
    typedef ug::ConvectionDiffusionPlugin::ConvectionDiffusionFV1<TDomain> base_type;
    typedef CplUserData<number, dim>                                       TNumberData;
    typedef CplUserData<MathVector<dim>, dim>                              TVectorData;
    typedef InverseLinker<dim>                                             TInverseLinker;
 
    StVenantElemDisc(const char* functions, const char* subsets)
        : base_type(functions, subsets),
          m_exValue(make_sp<DataExport<number, dim>>(new DataExport<number, dim>(functions))),
          m_exGrad(make_sp<DataExport<MathVector<dim>, dim>>(new DataExport<MathVector<dim>, dim>(functions)))
    {}
 
    // --- Getters and Setters ---
    SmartPtr<TNumberData>    get_width_data()  { return m_width; }
    SmartPtr<TInverseLinker> get_height_data() { return m_height; }
 
    void set_width_data(SmartPtr<TNumberData> b) { m_width = b; }
 
protected:
    SmartPtr<DataExport<number, dim>>             m_exValue;
    SmartPtr<DataExport<MathVector<dim>, dim>>    m_exGrad;
    SmartPtr<TNumberData>                         m_height;
    SmartPtr<TNumberData>                         m_width;
};
 
 
template<typename TDomain>
class StVenantElemDiscFactory {
public:
    static const int dim = TDomain::dim;
 
    // --- Type Definitions ---
    typedef IElemDisc<TDomain>                           TElemDisc;
    typedef StVenantElemDisc<TDomain>                    TStVenant;
    typedef ConstantSohlGradient<TDomain>                s_type;
    typedef EdgeOrientation<TDomain>                     TOrient;
    typedef ConvectionShapesFullUpwind<dim>              TUpwind;
 
    typedef CplUserData<number, dim>                     TNumberData;
    typedef CplUserData<MathVector<dim>, dim>            TVectorData;
    typedef ScaleAddLinker<MathVector<dim>, dim, number> TScaleAddLinkerVector;
    typedef ScaleAddLinker<number, dim, number>          TScaleAddLinkerNumber;
    typedef InverseLinker<dim>                           TInverseLinker;
 
    StVenantElemDiscFactory() {
        m_gravity = 9.81;
        m_diffusion = 1.0;
        m_functions = {"A", "v"}; // Cross-sectional Area (A) and Velocity (v)
    }
 
    virtual ~StVenantElemDiscFactory() = default;
 
    std::vector<SmartPtr<TStVenant>> create_1d_swe(const char* subsets) {
        UG_LOG("Creating 1D St. Venant element discretizations...\n");
 
        std::vector<SmartPtr<TStVenant>> elemDiscs;
 
        // 1. Initialize Base Discretizations
        SmartPtr<TStVenant> disc_A = make_sp<TStVenant>(new TStVenant(m_functions[0].c_str(), subsets));
        SmartPtr<TStVenant> disc_v = make_sp<TStVenant>(new TStVenant(m_functions[1].c_str(), subsets));
 
        auto A = disc_A->value();
        auto v = disc_v->value();
 
        // 2. Geometric Linkers
        // Water Depth: $h = A / b$
        auto h = make_sp(new TInverseLinker());
        h->divide(A, m_b);
 
        // 3. Stabilization (Upwinding)
        auto upwind = make_sp(new TUpwind());
        disc_A->set_upwind(upwind);
        disc_v->set_upwind(upwind);
 
        // 4. Advection/Velocity Vectors
        SmartPtr<TOrient> downStreamVec = make_sp(new TOrient(m_domain));
        auto boundaryDirection = make_sp(new StVenantWallBndDir<dim>(downStreamVec.template cast_dynamic<TVectorData>()));
        m_boundary_dir = boundaryDirection.template cast_dynamic<TNumberData>();
 
        // Convective velocity for the transport equations: $\vec{v}_{adv} = v \cdot \hat{d}$
        auto adv_velocity = make_sp(new TScaleAddLinkerVector());
        adv_velocity->add(v, downStreamVec);
 
        disc_A->set_velocity(adv_velocity);
        disc_v->set_velocity(adv_velocity);
 
        // 5. Flux/Pressure Terms for Velocity Equation
        // Absolute water surface elevation: $Z = h + z_{bed}$
        SmartPtr<s_type> bed_elevation = make_sp<s_type>(new s_type("Sink", 0.001, m_domain));
        m_s = bed_elevation;
 
 
 
        // Pressure gradient flux: $\vec{F} = g(h + z_{bed})\hat{d}$
        auto flux_inner = make_sp(new TScaleAddLinkerVector());
        flux_inner->add(h, downStreamVec);
        flux_inner->add(bed_elevation, downStreamVec);
 
        auto flux_v = make_sp(new TScaleAddLinkerVector());
        flux_v->add(m_gravity, flux_inner);
 
        m_v_flux_term = flux_v;
        disc_v->set_flux(flux_v);
 
        // 6. Diffusion
        disc_v->set_diffusion(m_diffusion);
 
        // 7. Reaction Rate / Friction Term
        auto roughness = make_sp(new ManningRoughness1d<dim>());
        roughness->set_h(h);
        roughness->set_v(v);
        roughness->set_manning_coefficient(m_frictionKs);
        roughness->set_hmin(1e-8);
        roughness->set_gravity(m_gravity);
 
        m_manningRoughness = roughness.template cast_dynamic<TNumberData>();
        disc_v->set_reaction_rate(roughness);
 
        // 8. Finalize
        elemDiscs.push_back(disc_A);
        elemDiscs.push_back(disc_v);
 
        UG_LOG("Successfully created St. Venant element discretizations.\n");
        return elemDiscs;
    }
 
    // --- Exporters ---
    SmartPtr<s_type>                export_sohle()             { return m_s; }
    SmartPtr<TScaleAddLinkerNumber> export_flux()              { return m_v_gsh; }
    SmartPtr<TScaleAddLinkerNumber> export_velocity()          { return m_v_05; }
    SmartPtr<TNumberData>           export_manning_roughness() { return m_manningRoughness; }
    SmartPtr<TVectorData>           export_v_flux_term()       { return m_v_flux_term; }
    SmartPtr<TScaleAddLinkerNumber> get_source_term()          { return m_total_source_term; }
    SmartPtr<TNumberData>           get_boundary_dir()         { return m_boundary_dir; }
    SmartPtr<TNumberData>           get_neumann_A()            { return m_A_neumann; }
    SmartPtr<TNumberData>           get_neumann_V()            { return m_v_neumann; }
 
    // --- Setters ---
    void set_domain(SmartPtr<TDomain> domain)         { m_domain = domain; }
    void set_gravity(number g)                        { m_gravity = g; }
    void set_b(SmartPtr<TNumberData> b)               { m_b = b; }
    void set_width(SmartPtr<TNumberData> b)           { m_b = b; } // Alias for set_b
    void set_ks(SmartPtr<TNumberData> ks)             { m_frictionKs = ks; }
    void set_diffusion(number nu)                     { m_diffusion = nu; }
    void set_h_threshold(number t)                    { m_threshold = t; }
 
    const char* get_disc_type() { return "fv1"; }
 
protected:
    std::vector<std::string>        m_functions;
    SmartPtr<TDomain>               m_domain;
 
    // Physical Parameters
    number                          m_gravity = 9.81;
    number                          m_diffusion = 0.5;
    number                          m_threshold = 1e-5;
 
    // User Data Fields
    SmartPtr<TNumberData>           m_frictionKs = make_sp<TNumberData>(new ConstUserNumber<dim>(0.1));
    SmartPtr<TNumberData>           m_b          = make_sp<TNumberData>(new ConstUserNumber<dim>(1.0));
    SmartPtr<s_type>                m_s;
 
    // Cached Linkers/Terms for export
    SmartPtr<TNumberData>           m_manningRoughness;
    SmartPtr<TVectorData>           m_v_flux_term;
    SmartPtr<TNumberData>           m_boundary_dir;
    SmartPtr<TNumberData>           m_A_neumann;
    SmartPtr<TNumberData>           m_v_neumann;
 
    SmartPtr<TScaleAddLinkerVector> m_v_flux     = make_sp<TScaleAddLinkerVector>(new TScaleAddLinkerVector());
    SmartPtr<TScaleAddLinkerNumber> m_v_gsh      = make_sp<TScaleAddLinkerNumber>(new TScaleAddLinkerNumber());
    SmartPtr<TScaleAddLinkerVector> m_A_velocity = make_sp<TScaleAddLinkerVector>(new TScaleAddLinkerVector());
    SmartPtr<TScaleAddLinkerVector> m_v_velocity = make_sp<TScaleAddLinkerVector>(new TScaleAddLinkerVector());
    SmartPtr<TScaleAddLinkerNumber> m_v_05       = make_sp<TScaleAddLinkerNumber>(new TScaleAddLinkerNumber());
    SmartPtr<TScaleAddLinkerNumber> m_total_source_term;
    SmartPtr<TUpwind>               m_upwind     = make_sp<TUpwind>(new TUpwind());
};
 
 
 
}
}
#endif //ST_VENANT_ELEM_DISC_H