add_river_sss.h

Go to the documentation of this file.

 
#ifndef UG4_ADD_RIVER_SSS_H
#define UG4_ADD_RIVER_SSS_H
 
#include "lib_disc/function_spaces/grid_function.h"
#include "lib_disc/domain_util.h"
#include "lib_disc/function_spaces/grid_function_global_user_data.h"
#include "../d3f_sss.h"
#include "lib_disc/function_spaces/integrate.h"
#include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
#include "lib_disc/spatial_disc/user_data/linker/scale_add_linker.h"
 
namespace ug {
namespace d3f {
 
template <typename TDomain>
class ExfiltrationFlux : public StdDataLinker<ExfiltrationFlux<TDomain>, number, TDomain::dim> {
public:
    static const int dim = TDomain::dim;
    typedef typename TDomain::grid_type TGrid;
    typedef typename TGrid::template EdgeAttachmentAccessor<ANumber> accessor_type;
 
    explicit ExfiltrationFlux(accessor_type& acc) : m_acc(acc) {}
 
    inline void evaluate(number& value, const MathVector<dim>& globIP, number time, int si) const {
        UG_THROW("ExfiltrationFlux: Point evaluation not supported. Use element-based evaluation.");
    }
 
    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
    {
        Edge* edgeElem = static_cast<Edge*>(elem);
        const number flux = m_acc[edgeElem];
        for (size_t ip = 0; ip < nip; ++ip) {
            vValue[ip] = flux;
        }
    }
 
    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
    {
        evaluate<refDim>(vValue, vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
        // Exchange flux is treated explicitly in the coupling; derivative w.r.t. surface state is zero.
    }
 
    bool requires_grid_fct() const { return true; }
 
private:
    accessor_type& m_acc;
};
 
template <typename TDomain, typename TAlgebra>
class RiverSSSManager {
public:
    static const int dim = TDomain::dim;
    typedef GridFunction<TDomain, TAlgebra>                            TGridFunction;
    typedef RiverSingularSourcesAndSinks<dim>                          TSSSManager;
    typedef typename TDomain::grid_type                                TGrid;
    typedef typename TGrid::template EdgeAttachmentAccessor<ANumber>   accessor_type;
    typedef typename TSSSManager::line_sss_type                        line_sss_type;
    typedef typename TGridFunction::element_type                       TElem;
 
    RiverSSSManager() = default;
 
    RiverSSSManager(SmartPtr<TGridFunction> riverU, SmartPtr<TGridFunction> ssU) 
        : m_spRiverU(riverU), m_spSubsurfaceU(ssU) 
    {
        m_spGlobEval = make_sp(new GlobalGridFunctionNumberData<TGridFunction>(ssU, "p"));
 
        TGrid& grid = *riverU->domain()->grid();
        grid.attach_to_edges(m_aExchangeFlux);
        m_aaExchangeFlux.access(grid, m_aExchangeFlux);
 
        m_spExfiltrationFlux = make_sp(new ExfiltrationFlux<TDomain>(m_aaExchangeFlux))
                                   .template cast_dynamic<CplUserData<number, dim>>();
    }
 
    // --- Configuration ---
    void set_bottom_permeability(number perm)      { m_perm = perm; }
    void set_bottom_thickness(number thickness)    { m_bottom_thickness = thickness; }
    void set_river_width(number width)             { m_width = width; }
    void set_density(SmartPtr<CplUserData<number,dim>> density) { m_spDensity = density; }
 
    SmartPtr<CplUserData<number, dim>> get_qex() { return m_spExfiltrationFlux; }
    void push_sss_manager(SmartPtr<TSSSManager> ssmngr) { m_vspSSSManager.push_back(ssmngr); }
 
    void init_line_sss() {
        UG_LOG("Initializing River SSS managers...\n");
        m_vspLineSSS.clear();
 
        auto river_grid = m_spRiverU->domain()->grid();
        const size_t numEdges = river_grid->template num<Edge>();
        m_vspLineSSS.assign(numEdges, SPNULL);
 
        auto vAcc = m_spRiverU->domain()->position_accessor();
        auto iter = m_spRiverU->template begin<Edge>();
        auto end  = m_spRiverU->template end<Edge>();
 
        size_t i = 0;
        for (; iter != end; ++iter, ++i) {
            Edge* edge = *iter;
            m_aaExchangeFlux[edge] = 0.0;
 
            MathVector<dim> x0 = vAcc[edge->vertex(0)];
            MathVector<dim> x1 = vAcc[edge->vertex(1)];
 
            // Evaluate subsurface pressure at the higher node (approximate phreatic surface)
            MathVector<dim> eval_point = (x0[dim - 1] > x1[dim - 1]) ? x0 : x1;
            
            number subsurface_pressure = 0.0;
            if (!m_spGlobEval->try_evaluate_global(subsurface_pressure, eval_point)) {
                // Processor doesn't own the 3D volume under this river segment
                continue; 
            }
 
            // Create SSS: Intensity is driven by the 1D water column height (h)
            auto spSSS = make_sp(new line_sss_type(x0, x1));
            spSSS->set_edge(edge);
            
            // We initialize with current height data
            std::vector<DoFIndex> ind;
            m_spRiverU->template dof_indices<Edge>(edge, 0, ind);
            number h = (x0[dim - 1] > x1[dim - 1]) ? DoFRef(*m_spRiverU, ind[0]) : DoFRef(*m_spRiverU, ind[1]);
 
            spSSS->set(m_gravity * m_density * h, 0.0, m_perm, m_width, m_bottom_thickness);
            m_vspLineSSS[i] = spSSS;
 
            for (auto& ssmngr : m_vspSSSManager) {
                ssmngr->add_line(spSSS);
            }
        }
    }
 
    void update_line_sss(number dt) {
        auto vAcc = m_spRiverU->domain()->position_accessor();
        auto iter = m_spRiverU->template begin<Edge>();
        auto end  = m_spRiverU->template end<Edge>();
 
        size_t i = 0;
        for (; iter != end; ++iter, ++i) {
            auto spSSS = m_vspLineSSS[i];
            if (!spSSS.valid()) continue;
 
            Edge* edge = *iter;
            std::vector<DoFIndex> ind;
            m_spRiverU->template dof_indices<Edge>(edge, 0, ind);
 
            MathVector<dim> x0 = vAcc[edge->vertex(0)];
            MathVector<dim> x1 = vAcc[edge->vertex(1)];
            number h = (x0[dim - 1] > x1[dim - 1]) ? DoFRef(*m_spRiverU, ind[0]) : DoFRef(*m_spRiverU, ind[1]);
 
            // Update SSS intensity: p_river = rho * g * h
            spSSS->set(m_gravity * m_density * h, 0.0, m_perm, m_width, m_bottom_thickness);
        }
    }
 
    void evaluate_exchange_flux(SmartPtr<TGridFunction> u) {
        const size_t numEdges = m_spRiverU->domain()->grid()->template num<Edge>();
        m_vflux.assign(numEdges, 0.0);
 
        // Reset local flux accumulation
        auto river_grid = m_spRiverU->domain()->grid();
        for (auto it = m_spRiverU->template begin<Edge>(); it != m_spRiverU->template end<Edge>(); ++it) {
            m_aaExchangeFlux[*it] = 0.0;
        }
 
        typedef typename TDomain::position_accessor_type t_pos_accessor;
        typedef typename RiverSingularSourcesAndSinks<dim>::template
            line_iterator<TElem, t_pos_accessor, DimFV1Geometry<dim>> t_lin_sss_iter;
 
        t_pos_accessor& aaPos = u->domain()->position_accessor();
        auto& grid = static_cast<Grid&>(*u->domain()->grid());
        DimFV1Geometry<dim> geo;
        std::vector<MathVector<dim>> vCorner;
 
        if (m_interface_elements.empty()) cache_interface_elements(u);
 
        // Iterate only over cached interface elements
        for (auto elem : m_interface_elements) {
            CollectCornerCoordinates(vCorner, *elem, aaPos, true);
            geo.update(elem, &vCorner[0], u->domain()->subset_handler().get());
 
            // Check each sub-control volume (SCV) of the 3D cell
            for (size_t ip = 0; ip < geo.num_scv(); ip++) {
                size_t co = geo.scv(ip).node_id();
                
                std::vector<DoFIndex> ind;
                u->dof_indices(elem, u->fct_id_by_name("p"), ind);
                const number p_sub = DoFRef(*u, ind[co]);
 
                // Find all river segments crossing this SCV
                for (t_lin_sss_iter line(m_vspSSSManager[0].get(), elem, grid, aaPos, geo, co); !line.is_over(); ++line) {
                    auto* line_sss = *line;
                    Edge* river_edge = line_sss->get_edge();
                    
                    const number segmentLength = VecDistance(line.seg_start(), line.seg_end());
                    const number p_riv        = line_sss->intensity();
                    const number matrix_perm  = line_sss->matrix_permeability();
                    const number rel_bot_perm = line_sss->relative_bottom_permeability();
 
                    // Only handle infiltration (river -> subsurface) if desired, 
                    // or both. Currently assumes exfiltration is driving the coupling.
                    if (p_sub < p_riv) continue;
 
                    // Reconstruct analytical leakage: Q [m^3/s] = L * (K/mu) * (p_sub - p_riv) * (W/B)
                    const number total_perm = matrix_perm * rel_bot_perm;
                    const number leakage_vol = segmentLength * (total_perm / m_viscosity) * 
                                              (p_sub - p_riv) * (m_width / m_bottom_thickness);
 
                    // Convert to target line flux [m^2/s]
                    const number edge_length = VecDistance(vAcc(river_edge->vertex(0)), vAcc(river_edge->vertex(1)));
                    m_aaExchangeFlux[river_edge] += leakage_vol / edge_length;
                }
            }
        }
 
        // --- MPI Synchronization ---
#ifdef UG_PARALLEL
        if (pcl::NumProcs() > 1) {
            // Flatten to array for reduction
            size_t idx = 0;
            for (auto it = m_spRiverU->template begin<Edge>(); it != m_spRiverU->template end<Edge>(); ++it, ++idx) {
                m_vflux[idx] = m_aaExchangeFlux[*it];
            }
            // Sum up partial fluxes from distributed 3D interface volumes
            pcl::ProcessCommunicator com;
            std::vector<number> global_fluxes(numEdges, 0.0);
            com.allreduce(&m_vflux[0], &global_fluxes[0], numEdges, pcl::DataTypeTraits<number>::get_data_type(), PCL_RO_SUM);
            // Write back global values to local 1D copies
            idx = 0;
            for (auto it = m_spRiverU->template begin<Edge>(); it != m_spRiverU->template end<Edge>(); ++it, ++idx) {
                m_aaExchangeFlux[*it] = global_fluxes[idx];
            }
        }
#endif
    }
 
private:
    void cache_interface_elements(SmartPtr<TGridFunction> u) {
        m_interface_elements.clear();
        typedef typename TDomain::position_accessor_type t_pos_accessor;
        typedef typename RiverSingularSourcesAndSinks<dim>::template
            line_iterator<TElem, t_pos_accessor, DimFV1Geometry<dim>> t_lin_sss_iter;
 
        t_pos_accessor& aaPos = u->domain()->position_accessor();
        auto& grid = static_cast<Grid&>(*u->domain()->grid());
        DimFV1Geometry<dim> geo;
        std::vector<MathVector<dim>> vCorner;
 
        for (auto iter = u->template begin<TElem>(); iter != u->template end<TElem>(); ++iter) {
            auto elem = *iter;
            CollectCornerCoordinates(vCorner, *elem, aaPos, true);
            geo.update(elem, &vCorner[0], u->domain()->subset_handler().get());
 
            bool found = false;
            for (size_t i = 0; i < geo.num_scv() && !found; i++) {
                t_lin_sss_iter line(m_vspSSSManager[0].get(), elem, grid, aaPos, geo, geo.scv(i).node_id());
                if (!line.is_over()) {
                    m_interface_elements.push_back(elem);
                    found = true;
                }
            }
        }
    }
 
    // Helper for vertex position access on the river grid
    MathVector<dim> vAcc(Vertex* v) const { return m_spRiverU->domain()->position_accessor()[v]; }
 
    SmartPtr<GlobalGridFunctionNumberData<TGridFunction>> m_spGlobEval;
    SmartPtr<TGridFunction>                               m_spRiverU;
    SmartPtr<TGridFunction>                               m_spSubsurfaceU;
    SmartPtr<CplUserData<number, dim>>                    m_spExfiltrationFlux;
    SmartPtr<CplUserData<number, dim>>                    m_spDensity;
 
    std::vector<SmartPtr<TSSSManager>>                    m_vspSSSManager;
    std::vector<SmartPtr<line_sss_type>>                  m_vspLineSSS;
    std::vector<number>                                   m_vflux;
    std::vector<TElem*>                                   m_interface_elements;
 
    number m_perm              = 1e-10;
    number m_width             = 1.0;
    number m_bottom_thickness  = 1.0;
    number m_gravity           = 9.81;
    number m_viscosity         = 1e-3;
    number m_density           = 998.2;
 
    ANumber         m_aExchangeFlux;
    accessor_type   m_aaExchangeFlux;
};
 
} // namespace d3f
} // namespace ug
 
#endif // UG4_ADD_RIVER_SSS_H