st_venant_submodule_impl.h

Go to the documentation of this file.

//
// Created by julian on 5/15/25.
//
 
#ifndef D3F_ST_VENANT_SUBMODULE_IMPL_H
#define D3F_ST_VENANT_SUBMODULE_IMPL_H
 
// --- Includes ---
#include <lib_algebra/operator/linear_solver/bicgstab.h>
#include <lib_algebra/operator/preconditioner/ilu.h>
#include <lib_disc/operator/non_linear_operator/newton_solver/newton.h>
#include <lib_disc/domain_util.h>
#include "st_venant_submodule.h"
#include <lib_disc/spatial_disc/domain_disc.h>
#include <lib_disc/spatial_disc/constraints/dirichlet_boundary/lagrange_dirichlet_boundary.h>
#include <lib_disc/spatial_disc/elem_disc/neumann_boundary/fv1/neumann_boundary_fv1.h>
#include "../../../Limex/time_disc/linear_implicit_timestep.h"
#include "../../../Limex/time_disc/simple_integrator.hpp"
#include "../../../Limex/time_disc/time_extrapolation.h"
#include "../../../Limex/time_disc/time_integrator.hpp"
#include "../../../Limex/time_disc/time_extrapolation.h"
#include "../../../Limex/time_disc/limex_integrator.hpp"
#include "lib_grid/file_io/file_io.h"
#include "lib_grid/file_io/file_io_ugx.h"
#include "../river_data_observer.h"
 
namespace ug {
namespace d3f {
 
template<typename TDomain, typename TAlgebra>
SmartPtr<LimexTimeIntegrator<TDomain, TAlgebra>> create_limex_solver(
    number tolerance,
    SmartPtr<DomainDiscretization<TDomain, TAlgebra>> domain_disc,
    SmartPtr<TDomain> domain,
    const char* cmp = "")
{
    // --- Type Definitions ---
    typedef typename TAlgebra::vector_type                   vector_type;
    typedef LimexTimeIntegrator<TDomain, TAlgebra>           limex_type;
    typedef NewtonSolver<TAlgebra>                           newton_solver_type;
    typedef BiCGStab<vector_type>                            linear_solver_type;
    typedef ILU<TAlgebra>                                    ilu_type;
    typedef StdConvCheck<vector_type>                        conv_check_type;
    typedef CompositeGridFunctionEstimator<TDomain, TAlgebra> gf_estimator_type;
    typedef IComponentSpace<GridFunction<TDomain,TAlgebra>>  composite_type;
    typedef L2ComponentSpace<GridFunction<TDomain,TAlgebra>> l2_type;
 
    // --- Core Solvers ---
    auto limex_integrator = make_sp(new limex_type(2));
    auto newton_solver    = make_sp(new newton_solver_type());
    auto linear_solver    = make_sp(new linear_solver_type());
    auto ilu              = make_sp(new ilu_type());
 
    // --- Convergence Checks ---
    auto linear_conv_check = make_sp(new conv_check_type(100, 1e-6, 1e-4));
    linear_solver->set_preconditioner(ilu);
    linear_solver->set_convergence_check(linear_conv_check);
 
    auto newton_conv_check = make_sp(new conv_check_type(1, 5e-12, 1e-12));
    newton_conv_check->set_supress_unsuccessful(true);
 
    newton_solver->set_linear_solver(linear_solver);
    newton_solver->set_convergence_check(newton_conv_check);
 
    // --- Error Estimators for Limex ---
    auto estimator = make_sp(new gf_estimator_type());
 
    // If a specific component is passed, estimate error only for that, otherwise do both A and v.
    if (cmp && cmp[0] != '\0') {
       auto l2_cmp = make_sp(new l2_type(cmp, 2));
       estimator->add(l2_cmp.template cast_dynamic<composite_type>());
    } else {
       auto l2_A = make_sp(new l2_type("A", 2));
       auto l2_v = make_sp(new l2_type("v", 2));
       estimator->add(l2_A.template cast_dynamic<composite_type>());
       estimator->add(l2_v.template cast_dynamic<composite_type>());
    }
 
    // --- Limex Configuration ---
    limex_integrator->add_stage(1, newton_solver, domain_disc);
    limex_integrator->add_stage(2, newton_solver, domain_disc);
 
    limex_integrator->set_tolerance(tolerance);
    limex_integrator->set_dt_min(1e-12);
    limex_integrator->set_dt_max(1e10);
    limex_integrator->add_error_estimator(estimator);
    limex_integrator->set_increase_factor(2.0);
    limex_integrator->select_cost_strategy(make_sp<LimexNonlinearCost>(new LimexNonlinearCost()));
 
    return limex_integrator;
}
 
template<typename TGridFunction, int dim>
void StVenantSubmodule<TGridFunction, dim>::run(number tStart, number tEnd, number dt) {
 
    // --- Type Definitions ---
    typedef typename TGridFunction::algebra_type                 algebra_type;
    typedef DomainDiscretization<domain_type, algebra_type>      domain_disc_type;
    typedef ScaleAddLinker<number, dim, number>                  scale_add_linker_type;
    typedef InverseLinker<dim>                                   inverse_linker_type;
    typedef NeumannBoundaryFV1<domain_type>                      neumann_type;
    typedef DirichletBoundary<domain_type, algebra_type>         dirichlet_type;
    typedef CplUserData<number, dim>                             TNumberData;
    typedef OutflowIntegrationObserver<domain_type, algebra_type> outflow_observer_type;
 
    // --- 1. Setup Discretization Factory ---
    auto factory = make_sp(new factory_type());
    factory->set_b(params.b);
    factory->set_gravity(params.gravity);
    factory->set_diffusion(params.diffusion);
    factory->set_domain(m_domain);
    factory->set_ks(params.ks);
 
    // Extract elemental discretizations for Cross-sectional Area (A) and Velocity (v)
    auto elem_discs = factory->create_1d_swe("River");
    m_elem_disc_A = elem_discs[0];
    m_elem_disc_v = elem_discs[1];
 
    auto A = m_elem_disc_A->value();
    auto v = m_elem_disc_v->value();
 
    // --- 2. Physics Linkers (Mathematical Expressions) ---
    // Water depth: h = A / b
    auto h = make_sp(new inverse_linker_type());
    h->divide(A, params.b);
 
    // Source terms for A: q_exfiltration + q_rainfall
    auto source_linker = make_sp(new scale_add_linker_type());
    source_linker->add(1.0, m_qex);
    source_linker->add(1.0, params.qr);
    m_elem_disc_A->set_source(source_linker);
 
    // --- 3. Boundary Conditions ---
    auto domain_disc = make_sp(new domain_disc_type(m_approx_space));
 
    auto neumann_A = make_sp(new neumann_type("A"));
    auto neumann_v = make_sp(new neumann_type("v"));
    auto dirichlet_A = make_sp(new dirichlet_type());
    auto dirichlet_v = make_sp(new dirichlet_type());
 
    // Flux out for Area: Flowrate Q = v * A
    auto flowrate = make_sp(new scale_add_linker_type());
    flowrate->add(v, A);
 
    // Flux out for Velocity: Includes gravity, water depth (h), riverbed elevation (sohle), and v^2
    auto vv = make_sp(new scale_add_linker_type());
    vv->add(v, v);
    auto outflow = make_sp(new scale_add_linker_type());
    outflow->add(params.gravity, h);
    outflow->add(params.gravity, factory->export_sohle());
    outflow->add(1.0, vv);
 
    // Directional linkers
    auto outflow_dir = factory->get_boundary_dir();
    auto wall_dir = make_sp(new scale_add_linker_type());
    wall_dir->add(-1.0, outflow_dir);
 
    auto out_outflow  = make_sp(new scale_add_linker_type());
    auto out_flowrate = make_sp(new scale_add_linker_type());
    out_outflow->add(outflow_dir, outflow);
    out_flowrate->add(outflow_dir, flowrate);
 
    // Assign BCs to respective subsets
    dirichlet_A->add(params.dirichlet_A, "A", "Level");
    dirichlet_v->add(0.0, "v", "Wall");
 
    neumann_A->add(out_flowrate.template cast_dynamic<TNumberData>(), "Sink", "River");
    neumann_v->add(out_outflow.template cast_dynamic<TNumberData>(), "Sink", "River");
    neumann_v->add(out_outflow.template cast_dynamic<TNumberData>(), "Level", "River");
 
    // Register all constraints and discretizations
    domain_disc->add(dirichlet_A.template cast_dynamic<IDomainConstraint<domain_type, algebra_type>>());
    domain_disc->add(dirichlet_v.template cast_dynamic<IDomainConstraint<domain_type, algebra_type>>());
    domain_disc->add(neumann_A.template cast_dynamic<IElemDisc<domain_type>>());
    domain_disc->add(neumann_v.template cast_dynamic<IElemDisc<domain_type>>());
    domain_disc->add(m_elem_disc_A.template cast_dynamic<IElemDisc<domain_type>>());
    domain_disc->add(m_elem_disc_v.template cast_dynamic<IElemDisc<domain_type>>());
 
    // --- 4. Observers and VTK Output ---
    auto limex_integrator = create_limex_solver(m_tol, domain_disc, this->m_u->domain());
    auto downStreamVec = make_sp(new EdgeOrientation<domain_type>(m_domain));
 
    auto vtk = make_sp<VTKOutput<dim>>(new VTKOutput<dim>());
    vtk->select("v", "v");
    vtk->select("A", "A");
    vtk->select(factory->export_sohle().template cast_dynamic<UserData<number, dim>>(), "sohle");
    vtk->select(flowrate.template cast_dynamic<UserData<number, dim>>(), "flowrate");
    vtk->select(downStreamVec.template cast_dynamic<UserData<MathVector<dim>, dim>>(), "EdgeOrientation");
    vtk->select(m_qex.template cast_dynamic<UserData<number, dim>>(), "qex");
    auto voutflow = make_sp(new ScaleAddLinker<MathVector<dim>, dim, number>());
    voutflow->add(v, downStreamVec);
    auto outflow_observer = make_sp(new outflow_observer_type(voutflow, A, "Sink", "River"));
    limex_integrator->attach_observer(outflow_observer);
 
    // --- 5. Execution ---
#ifdef UG_PARALLEL
    pcl::ProcessCommunicator com;
    com.barrier();
    UG_LOG_ALL_PROCS("starting limex\n");
#endif
 
    UG_LOG_ALL_PROCS("integrating from " << tStart << " to " << tEnd << "\n");
 
    if (!m_domain->empty() || true) {
        limex_integrator->apply(this->m_u, tEnd, this->m_u, tStart);
    }
 
    #ifdef UG_PARALLEL
    com.barrier();
    if (pcl::ProcRank() == 0 and m_step > 0) {
    #endif
 
        typedef typename TGridFunction::template dim_traits<1>::const_iterator const_iterator;
        typedef typename domain_traits<1>::grid_base_object grid_base_object;
        const_iterator iter = this->m_u->template begin<grid_base_object>();
        const_iterator iterEnd = this->m_u->template end<grid_base_object>();
        typedef UserDataIntegrand<number,TGridFunction> integrand_type;
        int quad_order = 1;
 
        integrand_type A_integrand = integrand_type(m_elem_disc_A->value(), this->m_u.get(), tEnd);
        number value = 1000 * Integrate<dim, 1, const_iterator>(iter, iterEnd,
                                                         this->m_u->domain()->position_accessor(),
                                                         A_integrand,
                                                         quad_order, "");
        vtk->print("results/rivers", *this->m_u, m_step, tEnd);
 
        std::ofstream myfile;
        myfile.open("results/river_outflow.txt", std::ios::out | std::ios::app);
        myfile << tEnd << "\t" << outflow_observer->value() << "\t"<< value << "\n";
        UG_LOGN("Total River outflow for dt = " <<tEnd - tStart << " is " << outflow_observer->value() << " kg\n" );
        myfile.close();
#ifdef UG_PARALLEL
 
    }
#endif
 
    m_step++;
}
 
template<typename TGridFunction, int dim>
void StVenantSubmodule<TGridFunction, dim>::init(const char* filename) {
 
    const char* m_filename  = "results/river_outflow.txt";
    #ifdef UG_PARALLEL
    if (pcl::ProcRank() == 0) 
    #endif
    {
        std::ofstream myfile;
        myfile.open(m_filename,  std::ios::out | std::ios::trunc);
        myfile << "time\taccumulated outflow since last measurement\tstorage\n";
        myfile.close();
    }
 
    // --- 1. Setup Domain & Grid Load ---
    SmartPtr<typename TGridFunction::domain_type> dom = make_sp(new typename TGridFunction::domain_type());
    dom->update_domain_info();
 
    std::vector<std::string> additionalSHNames = dom->additional_subset_handler_names();
    std::vector<SmartPtr<ISubsetHandler>> ash(additionalSHNames.size());
    for(size_t i = 0; i < additionalSHNames.size(); ++i) {
       ash[i] = dom->additional_subset_handler(additionalSHNames[i]);
    }
 
    dom->grid()->message_hub()->post_message(GridMessage_Creation(GMCT_CREATION_STARTS, -2));
 
    SPProjectionHandler ph = make_sp(new ProjectionHandler(dom->geometry3d(), dom->subset_handler()));
    size_t num_ph = 0;
 
    // Load mesh from UGX file
    LoadGridFromUGX(*(dom->grid()), ph, num_ph, *(dom->subset_handler()), additionalSHNames, ash, filename, dom->position_attachment());
 
    dom->update_subset_infos(-2);
    dom->grid()->message_hub()->post_message(GridMessage_Creation(GMCT_CREATION_STOPS));
    dom->update_domain_info();
 
    m_domain = dom;
 
    // --- 2. Setup Approximation Space ---
    m_approx_space = make_sp(new approx_space_type(m_domain));
 
    // Initialize Lagrange P1 functions for Cross-Sectional Area (A) and Velocity (v)
    m_approx_space->add("A", "Lagrange", 1);
    m_approx_space->add("v", "Lagrange", 1);
    m_approx_space->init_levels();
    m_approx_space->init_top_surface();
 
    // --- 3. Initial Interpolations ---
    this->m_u = make_sp(new TGridFunction(m_approx_space));
    Interpolate(params.dirichlet_A, this->m_u, "A");
    Interpolate(0.0, this->m_u, "v");
 
    // --- 4. Initialize Discretizations ---
    auto factory = make_sp(new factory_type());
    factory->set_domain(m_domain);
    factory->set_b(params.b);
    factory->set_gravity(params.gravity);
    factory->set_diffusion(params.diffusion);
    factory->set_ks(params.ks);
 
    auto elem_discs = factory->create_1d_swe("River");
    m_elem_disc_A = elem_discs[0];
    m_elem_disc_v = elem_discs[1];
}
 
} // namespace d3f
} // namespace ug
 
#endif //D3F_ST_VENANT_SUBMODULE_IMPL_H