fs_equilibirum_finished_condition.hpp

Go to the documentation of this file.

/* 
 * SPDX-FileCopyrightText: 2025 Gesellschaft fuer Anlagen- und Reaktorsicherheit gGmbH
 * SPDX-License-Identifier: EUPL-1.2
 * SPDX-FileContributor: Tim Schoen
 * SPDX-FileContributor: Goethe Universität Frankfurt
 * SPDX-FileType: SOURCE
 * 
 * This file is part of d3f++.
 * d3f++ is an extension for UG4. Licensing information and citation requirements of UG4 are provided in LICENSES/UG4-LGPL_2.1
 */
 
 
 
#ifndef __H__UG__PLUGINS__D3F__FS_EQUIBILIRIUM_FINISHED_CONDITION
#define __H__UG__PLUGINS__D3F__FS_EQUIBILIRIUM_FINISHED_CONDITION
 
 
#include "level_set_pos.h"
#include <lib_disc/time_disc/time_integrator_observers/time_integrator_observer_interface.h>
#include <lib_disc/time_disc/finished_conditions.hpp>
#include <fstream>
#include <cmath>
 
#ifdef UG_PARALLEL
#include "pcl/pcl_base.h"
#include "pcl/pcl_process_communicator.h"
#endif
 
namespace ug {
 
 
template<class TDomain, class TAlgebra>
class FSInEquilibriumFinishedCondition : public ITimeIntegratorObserver<TDomain, TAlgebra>, public IFinishedCondition
{
    private:
 
 
        static const int dim = TDomain::dim;
        struct t_pnt_data
        {
            MathVector<dim-1> xy; 
            number z; 
            
            bool valid; 
            
            t_pnt_data (MathVector<dim-1> & the_xy)
            : xy (the_xy), z (0), valid (false) {};
        };
 
 
    public:
 
      typedef ug::GridFunction<TDomain, TAlgebra> grid_function_type;
 
        FSInEquilibriumFinishedCondition(SmartPtr<grid_function_type> lsf, number equilibrium_factor) : m_equilibrium_factor(equilibrium_factor), m_height_measurer(lsf)
        {
          m_max_change = 0;
            m_last_change = 0;
        }
 
        void set_equilibrium_factor(number equilibrium_factor)
        {
            m_equilibrium_factor = equilibrium_factor;
        }
 
        bool check_finished(number time, int step)
        {        
            bool finished = false;
            
            #ifdef UG_PARALLEL
            if(pcl::ProcRank() == 0)
            {
            #endif
 
            if(m_last_measurement.size() == 0)             
            {
                // no measurement has happenend yet
                std::cout << "FSInEquilibriumFinishedCondition: no measurement has happenend yet" << std::endl;
                finished = false;
            }
            else if(m_max_change == 0)
            {
                std::cout << "FSInEquilibriumFinishedCondition: max_change is 0" << std::endl;
                finished = false;
            }
            else
            {                
                std::cout << "FSInEquilibriumFinishedCondition: factor is " << m_last_change << " / " << m_max_change << " = " << m_last_change/m_max_change << std::endl;
                finished = m_last_change/m_max_change < m_equilibrium_factor;
            }            
 
 
            #ifdef UG_PARALLEL
            }
      pcl::ProcessCommunicator proc_comm;
      proc_comm.broadcast(finished, 0);
            
            #endif
 
                
            return finished;
        }
 
        void add_measurement_point(std::vector<number> point)
        {        
            if(point.size() != TDomain::dim-1)
            {
                std::cout << "FSInEquilibriumFinishedCondition: free surface measurement point does not have correct dimension " << TDomain::dim-1 << ", ignoring!" << std::endl;
                return;
            }
 
            MathVector<dim-1> p;
 
            for(int i = 0; i < TDomain::dim-1; i++)
            {
                p[i] = point[i];
            }
                           
            m_measurement_points.push_back(t_pnt_data(p));
            
        }
 
      virtual bool step_process(SmartPtr<grid_function_type> u, int step, number time, number dt)
        {
 
            m_height_measurer.reinit();
 
            #ifndef UG_PARALLEL
 
            for (size_t i = 0; i < m_measurement_points.size (); i++)
            {
                t_pnt_data & pnt = m_measurement_points[i];
                pnt.valid = m_height_measurer.get_height_at (pnt.xy, pnt.z);
            }
  
            #else // i.e. ifdef UG_PARALLEL
                
            {
 
                pcl::ProcessCommunicator proc_comm;
                std::vector<number> loc_z (m_measurement_points.size ()), red_z (m_measurement_points.size ());
                
                for (size_t i = 0; i < m_measurement_points.size (); i++)
                    if (! m_height_measurer.get_height_at (m_measurement_points[i].xy, loc_z[i]))
                        loc_z[i] = - std::numeric_limits<number>::max (); // the position is not covered by the free surface in this process
                
                proc_comm.allreduce (loc_z, red_z, PCL_RO_MAX);
                for (size_t i = 0; i < m_measurement_points.size (); i++)
                {
                    t_pnt_data & pnt = m_measurement_points[i];
                    number z_val = red_z[i];
                    if (z_val != - std::numeric_limits<number>::max ()) // i.e. is not set in any process
                    {
                        pnt.z = z_val;
                        pnt.valid = true;
                    }
                    else
                        pnt.valid = false;
                }
            }
                
            if(pcl::ProcRank() != 0)
                return true;
            #endif
 
 
            //##############################################################
            std::vector<number> this_measurement;
 
            // do the needed measurements
            for (size_t i = 0; i < m_measurement_points.size(); i++)
            {
                number height = m_measurement_points[i].z;
 
                if (!m_measurement_points[i].valid)
                    height = 0;
 
                this_measurement.push_back(height);
            }
 
            if (m_last_measurement.size() == 0)
            {
                m_last_measurement = this_measurement;
                return true;
            }
 
            if (m_last_measurement.size() != this_measurement.size())
            {
                throw std::logic_error("Number of measurement points changed during runtime!");
            }         
 
 
            // calculate the 2-norm of the difference vector
            m_last_change = 0;
 
            for (size_t i = 0; i < m_last_measurement.size(); i++)
            {              
                if (m_measurement_points[i].valid)
                {
                    double difference = m_last_measurement[i]-this_measurement[i];
                    m_last_change += difference*difference;                    
                }
                else
                {
                    std::cout << "EquilibriumFinishedTester: No free surface at " << m_measurement_points[i].xy[0] << std::endl;
                }                
            }
 
            m_last_change = sqrt(m_last_change);
            m_last_change /= dt;
               
            if(m_last_change > m_max_change)
            {
                m_max_change = m_last_change;
            }
 
            m_last_measurement = this_measurement;
 
            return true;
 
        }
 
        void read_measurement_points(const char * filename)
        {
 
            std::string full_file_name = FindFileInStandardPaths (filename);
            UG_COND_THROW (full_file_name.empty (), "FSInEquilibriumFinishedCondition: Coulnd't locate file '" << filename << "'.");
            
            std::ifstream input (full_file_name.c_str (), std::ifstream::in);
            UG_COND_THROW (input.fail (), "FSInEquilibriumFinishedCondition: Cannot open data file '" << full_file_name << "' for input!");
 
            std::string input_line;
            std::istringstream line_stream;
            
            line_stream.exceptions (std::istream::failbit | std::istream::badbit);
            
            while (! input.eof ())
            {
                std::vector<number> xy(dim-1);
                
                if (input.fail ())
                    UG_THROW ("FSInEquilibriumFinishedCondition: Could not load the points from the file!");
                std::getline (input, input_line);
                if (input_line.length () == 0)
                    continue;
                try
                {
                    line_stream.str (input_line);
                    line_stream.clear ();
                    for (size_t i = 0; i < dim-1; i++)
                        line_stream >> xy[i];
                }
                catch (std::istream::failure & e)
                {
                    UG_THROW ("FSInEquilibriumFinishedCondition: Failed to parse line '" << input_line << "' in the input file.");
                };
            
                add_measurement_point (xy);
            }
        }
 
    private:
        int m_number_of_measurement_points;
        number m_equilibrium_factor;
        number m_max_change;
        number m_last_change;
        std::vector<t_pnt_data> m_measurement_points;
        std::vector<number> m_last_measurement;
        ug::d3f::LSPositionZ<grid_function_type> m_height_measurer;        
 
};
 
}
 
#endif