/******************************************************************************
 * @internal
 * @file        rCFD_balance.h
 * @copyright   (C) 2023
 *              Particulate Flow Modelling,
 *              Johannes Kepler University, Linz, Austria
 * @website     www.particulate-flow.at
 * @authors     Stefan Pirker
 *****************************************************************************/

#ifndef RCFD_BALANCE
#define RCFD_BALANCE

#include "rCFD_types.h"
#include "rCFD_macros.h"
#include "rCFD_topo.h"

    void balance_update_mass(rCFD_type* rcfd, const int i_layer, const int i_phase)
    {
#if RP_NODE

#if 1   // local vars
        int     i_data, i_cell, i_frame;
#endif

        // init. Balance.mass_integral
        for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

            rcfd->Balance[i_phase][i_data].mass_integral = 0.0;
        }

        i_frame = rcfd->Rec->global_frame;

        // accumulate Balance.mass_integral
        for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

#ifdef RCFD_LIB_ISLANDS
            i_frame = rcfd->Rec->global_frame_of_island[_C(i_layer).island_id[i_cell]];
#endif

            for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

                switch (rcfd->Data_Dict[i_phase][i_data].type){

                    case RCFD_DATA_TEMPERATURE:
                    {
                        rcfd->Balance[i_phase][i_data].mass_integral +=   _C(i_layer).data[i_phase][i_cell][i_data]
                                                                        * _C(i_layer).volume[i_cell]
                                                                        * _C(i_layer).vof[i_frame][i_cell][i_phase]
                                                                        * rcfd->Phase_Dict[i_phase].density
                                                                        * rcfd->Phase_Dict[i_phase].heat_capacity;
                        break;
                    }
                    case RCFD_DATA_CONCENTRATION:
                    {
                        rcfd->Balance[i_phase][i_data].mass_integral +=   _C(i_layer).data[i_phase][i_cell][i_data]
                                                                        * _C(i_layer).volume[i_cell]
                                                                        * _C(i_layer).vof[i_frame][i_cell][i_phase]
                                                                        * rcfd->Phase_Dict[i_phase].density;
                        break;
                    }
                    case RCFD_DATA_GENERIC:
                    {
                        rcfd->Balance[i_phase][i_data].mass_integral += _C(i_layer).data[i_phase][i_cell][i_data];
                        break;
                    }
                    default:
                        break;
                }
            }
        }

        // calc. Balance.mass_integral_global
        for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

            rcfd->Balance[i_phase][i_data].mass_integral_global = PRF_GRSUM1(rcfd->Balance[i_phase][i_data].mass_integral);
        }

#endif
    }

    void balance_init_source(rCFD_type* rcfd, const int i_layer, const int i_phase)
    {
#if RP_NODE

#if 1   // local vars
        int     i_data;
#endif

        // reset Balance.mass_source
        for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

            rcfd->Balance[i_phase][i_data].mass_source = 0.0;
            rcfd->Balance[i_phase][i_data].mass_source_global = 0.0;
        }

#endif
    }

    void balance_update_source(rCFD_type* rcfd, const int i_layer, const int i_phase)
    {
#if RP_NODE

#if 1   // local vars
        int     i_data;
#endif

        // add sources to Balance.mass_integral_target
        for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

            rcfd->Balance[i_phase][i_data].mass_integral_target += rcfd->Balance[i_phase][i_data].mass_source;
            rcfd->Balance[i_phase][i_data].mass_source_global = PRF_GRSUM1(rcfd->Balance[i_phase][i_data].mass_source);
            rcfd->Balance[i_phase][i_data].mass_integral_target_global += rcfd->Balance[i_phase][i_data].mass_source_global;
        }
#endif
    }

    void balance_init_cross_partition_fluxes(rCFD_type* rcfd, const int i_layer, const int i_phase)
    {
#if RP_NODE

#if 1   // local vars
        int     i_data, i_node, i_node2;
#endif

        // init Balance.node2node_flux
        for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

            if(rcfd->Balance_Dict[i_phase][i_data].type == RCFD_BALANCE_LOCAL){

                for(i_node = 0; i_node < compute_node_count; i_node++){

                    for(i_node2 = 0; i_node2 < compute_node_count; i_node2++){

                        rcfd->Balance[i_phase][i_data].node2node_flux[i_node][i_node2] = 0.0;
                        rcfd->Balance[i_phase][i_data].node2node_data_flux[i_node][i_node2] = 0.0;
                    }
                }
            }
        }
#endif
    }

    void balance_update_cross_partition_fluxes(rCFD_type* rcfd, const int i_layer, const int i_phase)
    {
#if RP_NODE
#endif
    }

    void balance_update_error(rCFD_type* rcfd, const int i_layer, const int i_phase)
    {
#if RP_NODE

#if 1   // local vars
        int     i_data;
#endif

        // calc. Balance.mass_error
        for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

            rcfd->Balance[i_phase][i_data].mass_error =   rcfd->Balance[i_phase][i_data].mass_integral_target
                                                        - rcfd->Balance[i_phase][i_data].mass_integral;
            rcfd->Balance[i_phase][i_data].mass_error_global = PRF_GRSUM1(rcfd->Balance[i_phase][i_data].mass_error);
        }
#endif
    }

    void balance_update(rCFD_type* rcfd, const int i_phase)
    {
#if RP_NODE
        int i_layer = rcfd->Solver->current_layer;

        balance_update_mass(rcfd, i_layer, i_phase);
        balance_update_source(rcfd, i_layer, i_phase);
        balance_update_cross_partition_fluxes(rcfd, i_layer, i_phase);
        balance_update_error(rcfd, i_layer, i_phase);
#endif
    }

    void balance_correct(rCFD_type* rcfd, const int i_phase)
    {
#if RP_NODE

    // TODO: Avoid code duplications for different data types!!

#if 1   // local vars
        int     i_layer, i_data, i_while, i_frame, i_cell, i_face;
        int     balance_error_exists, c0, c1, c, i_frame_c0, i_frame_c1;
        double  available_exchange, vol_flip, exchange_ratio;
#endif

        i_layer = rcfd->Solver->current_layer;

        if((rcfd->Solver->global_run_counter % rcfd->Solver_Dict->balance_correction_update) == 0){

            for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

                switch(rcfd->Data_Dict[i_phase][i_data].type){

                    case RCFD_DATA_GENERIC:
                    {
                        i_while = 0;

                        balance_error_exists = 1;

                        while((balance_error_exists) && (i_while < rcfd->Balance_Dict[i_phase][i_data].max_correction_loops)){

                            for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

                                _C(i_layer).data_swap[i_phase][i_cell][i_data] = 0.0;
                            }

                            available_exchange = 0.0;

                            for(i_face = 0; i_face < _FACE_DICT(i_layer).number_of_int_faces; i_face++){

                                c0 = _F(i_layer).c0[i_face];
                                c1 = _F(i_layer).c1[i_face];

                                available_exchange +=   0.5
                                                      * fabs(  _C(i_layer).data[i_phase][c1][i_data]
                                                             - _C(i_layer).data[i_phase][c0][i_data])
                                                      * rcfd->Solver_Dict->face_swap_max_per_loop;
                            }

                            if(rcfd->Balance_Dict[i_phase][i_data].type == RCFD_BALANCE_GLOBAL){

                                available_exchange = PRF_GRSUM1(available_exchange);
                            }

                            if(available_exchange > fabs(rcfd->Balance[i_phase][i_data].mass_error_global)){

                                if(available_exchange > 0.0){

                                    exchange_ratio = fabs(rcfd->Balance[i_phase][i_data].mass_error_global) / available_exchange;
                                }
                                else{

                                    exchange_ratio = 0.0;
                                }
                            }
                            else{

                                exchange_ratio = 1.0;
                            }

                            for(i_face = 0; i_face < _FACE_DICT(i_layer).number_of_int_faces; i_face++){

                                c0 = _F(i_layer).c0[i_face];
                                c1 = _F(i_layer).c1[i_face];

                                if(_C(i_layer).data[i_phase][c0][i_data] > _C(i_layer).data[i_phase][c1][i_data]){

                                    c = c0, c1 = c0, c0 = c;
                                }

                                if(_C(i_layer).data[i_phase][c0][i_data] < _C(i_layer).data[i_phase][c1][i_data]){

                                    if(rcfd->Balance[i_phase][i_data].mass_error_global > 0.0){

                                        _C(i_layer).data_swap[i_phase][c0][i_data] += 0.5
                                                                                    * exchange_ratio
                                                                                    * (_C(i_layer).data[i_phase][c1][i_data] - _C(i_layer).data[i_phase][c0][i_data])
                                                                                    * rcfd->Solver_Dict->face_swap_max_per_loop;
                                    }

                                    if(rcfd->Balance[i_phase][i_data].mass_error_global < 0.0){

                                        _C(i_layer).data_swap[i_phase][c1][i_data] -= 0.5
                                                                                    * exchange_ratio
                                                                                    * (_C(i_layer).data[i_phase][c1][i_data] - _C(i_layer).data[i_phase][c0][i_data])
                                                                                    * rcfd->Solver_Dict->face_swap_max_per_loop;
                                    }
                                }
                            }

                            rcfd->Balance[i_phase][i_data].mass_integral = 0.0;

                            for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

                                _C(i_layer).data[i_phase][i_cell][i_data] += _C(i_layer).data_swap[i_phase][i_cell][i_data];
                                rcfd->Balance[i_phase][i_data].mass_integral += _C(i_layer).data[i_phase][i_cell][i_data];
                            }

                            rcfd->Balance[i_phase][i_data].mass_integral_global = PRF_GRSUM1(rcfd->Balance[i_phase][i_data].mass_integral);
                            rcfd->Balance[i_phase][i_data].mass_error_global =    rcfd->Balance[i_phase][i_data].mass_integral_target_global
                                                                                - rcfd->Balance[i_phase][i_data].mass_integral_global;

                            if(  fabs(rcfd->Balance[i_phase][i_data].mass_error_global)
                               < fabs(  rcfd->Balance_Dict[i_phase][i_data].accuracy_level
                                      * rcfd->Balance[i_phase][i_data].mass_integral_target_global)){

                                balance_error_exists = 0;
                            }

                            i_while++;
                        }

                        if(i_while == rcfd->Balance_Dict[i_phase][i_data].max_correction_loops){

                            report_warning0(rcfd, RCFD_REPORT_CATEGORY_RUN_BALANCE,
                                            "Balance for i_phase %d i_data %d beyond accuracy level %f",
                                            i_phase, i_data, rcfd->Balance_Dict[i_phase][i_data].accuracy_level);
                        }

                        break;
                    }
                    case RCFD_DATA_CONCENTRATION:
                    {
                        i_while = 0;

                        balance_error_exists = 1;

                        i_frame_c0 = rcfd->Rec->global_frame;
                        i_frame_c1 = rcfd->Rec->global_frame;

                        while((balance_error_exists) && (i_while < rcfd->Balance_Dict[i_phase][i_data].max_correction_loops)){

                            // init. _C(i_layer).data_swap
                            for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

                                _C(i_layer).data_swap[i_phase][i_cell][i_data] = 0.0;
                            }

                            available_exchange = 0.0;

                            for(i_face = 0; i_face < _FACE_DICT(i_layer).number_of_int_faces; i_face++){

                                c0 = _F(i_layer).c0[i_face];
                                c1 = _F(i_layer).c1[i_face];

                                if(_C(i_layer).data[i_phase][c0][i_data] > _C(i_layer).data[i_phase][c1][i_data]){

                                    c = c0; c0 = c1; c1 = c;
                                }

                                if(_C(i_layer).data[i_phase][c1][i_data] > _C(i_layer).data[i_phase][c0][i_data]){

#ifdef RCFD_LIB_ISLANDS
                                    i_frame_c0 = rcfd->Rec->global_frame_of_island[_C(i_layer).island_id[c0]];
                                    i_frame_c1 = rcfd->Rec->global_frame_of_island[_C(i_layer).island_id[c1]];
#endif

                                    if(   (_C(i_layer).volume[c0] * _C(i_layer).vof[i_frame_c0][c0][i_phase])
                                        < (_C(i_layer).volume[c1] * _C(i_layer).vof[i_frame_c1][c1][i_phase])){

                                        vol_flip = _C(i_layer).volume[c0] * _C(i_layer).vof[i_frame_c0][c0][i_phase];
                                    }
                                    else{

                                        vol_flip = _C(i_layer).volume[c1] * _C(i_layer).vof[i_frame_c1][c1][i_phase];
                                    }

                                    available_exchange += 0.5
                                                        * (_C(i_layer).data[i_phase][c1][i_data] - _C(i_layer).data[i_phase][c0][i_data])
                                                        * vol_flip * rcfd->Phase_Dict[i_phase].density
                                                        * rcfd->Solver_Dict->face_swap_max_per_loop;
                                }
                            }

                            if(rcfd->Balance_Dict[i_phase][i_data].type == RCFD_BALANCE_GLOBAL){

                                available_exchange = PRF_GRSUM1(available_exchange);

                                if(available_exchange > fabs(rcfd->Balance[i_phase][i_data].mass_error_global)){

                                    if(available_exchange > 0.0){

                                        exchange_ratio = fabs(rcfd->Balance[i_phase][i_data].mass_error_global)
                                                        / available_exchange;
                                    }
                                    else{

                                        exchange_ratio = 0.0;
                                    }
                                }
                                else{

                                    exchange_ratio = 1.0;
                                }
                            }
                            else{

                                exchange_ratio = 0.0;
                            }

                            for(i_face = 0; i_face < _FACE_DICT(i_layer).number_of_int_faces; i_face++){

                                c0 = _F(i_layer).c0[i_face];
                                c1 = _F(i_layer).c1[i_face];

                                if(_C(i_layer).data[i_phase][c0][i_data] > _C(i_layer).data[i_phase][c1][i_data]){

                                    c = c0, c1 = c0, c0 = c;
                                }

                                if(_C(i_layer).data[i_phase][c0][i_data] < _C(i_layer).data[i_phase][c1][i_data]){

#ifdef RCFD_LIB_ISLANDS
                                    i_frame_c0 = rcfd->Rec->global_frame_of_island[_C(i_layer).island_id[c0]];
                                    i_frame_c1 = rcfd->Rec->global_frame_of_island[_C(i_layer).island_id[c1]];
#endif

                                    if((_C(i_layer).volume[c0] * _C(i_layer).vof[i_frame_c0][c0][i_phase]) < (_C(i_layer).volume[c1] * _C(i_layer).vof[i_frame_c1][c1][i_phase])){

                                        vol_flip = _C(i_layer).volume[c0] * _C(i_layer).vof[i_frame_c0][c0][i_phase];
                                    }
                                    else{
                                        vol_flip = _C(i_layer).volume[c1] * _C(i_layer).vof[i_frame_c1][c1][i_phase];
                                    }

                                    if(rcfd->Balance[i_phase][i_data].mass_error_global > 0.0){

                                        _C(i_layer).data_swap[i_phase][c0][i_data] += exchange_ratio
                                                                                    * (_C(i_layer).data[i_phase][c1][i_data] - _C(i_layer).data[i_phase][c0][i_data])
                                                                                    * vol_flip
                                                                                    * rcfd->Phase_Dict[i_phase].density
                                                                                    * rcfd->Solver_Dict->face_swap_max_per_loop;
                                    }

                                    if(rcfd->Balance[i_phase][i_data].mass_error_global < 0.0){

                                        _C(i_layer).data_swap[i_phase][c1][i_data] -= exchange_ratio
                                                                                    * (_C(i_layer).data[i_phase][c1][i_data] - _C(i_layer).data[i_phase][c0][i_data])
                                                                                    * vol_flip
                                                                                    * rcfd->Phase_Dict[i_phase].density
                                                                                    * rcfd->Solver_Dict->face_swap_max_per_loop;
                                    }
                                }
                            }

                            rcfd->Balance[i_phase][i_data].mass_integral = 0.0;

                            i_frame = rcfd->Rec->global_frame;

                            for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

#ifdef RCFD_LIB_ISLANDS
                                i_frame = rcfd->Rec->global_frame_of_island[_C(i_layer).island_id[i_cell]];
#endif

                                if((_C(i_layer).volume[i_cell] * _C(i_layer).vof[i_frame][i_cell][i_phase] * rcfd->Phase_Dict[i_phase].density) > 0.0){

                                    _C(i_layer).data[i_phase][i_cell][i_data] = (
                                                                                    _C(i_layer).data[i_phase][i_cell][i_data]
                                                                                  * _C(i_layer).volume[i_cell]
                                                                                  * _C(i_layer).vof[i_frame][i_cell][i_phase]
                                                                                  * rcfd->Phase_Dict[i_phase].density
                                                                                  + _C(i_layer).data_swap[i_phase][i_cell][i_data]
                                                                                )
                                                                                /
                                                                                (
                                                                                    _C(i_layer).volume[i_cell]
                                                                                  * _C(i_layer).vof[i_frame][i_cell][i_phase]
                                                                                  * rcfd->Phase_Dict[i_phase].density
                                                                                );
                                }

                                rcfd->Balance[i_phase][i_data].mass_integral +=   _C(i_layer).data[i_phase][i_cell][i_data]
                                                                                * _C(i_layer).volume[i_cell]
                                                                                * _C(i_layer).vof[i_frame][i_cell][i_phase]
                                                                                * rcfd->Phase_Dict[i_phase].density;
                            }

                            rcfd->Balance[i_phase][i_data].mass_integral_global = PRF_GRSUM1(rcfd->Balance[i_phase][i_data].mass_integral);
                            rcfd->Balance[i_phase][i_data].mass_error_global =    rcfd->Balance[i_phase][i_data].mass_integral_target_global
                                                                                - rcfd->Balance[i_phase][i_data].mass_integral_global;

                            if(fabs(rcfd->Balance[i_phase][i_data].mass_error_global)
                                < (  rcfd->Balance_Dict[i_phase][i_data].accuracy_level
                                   * rcfd->Balance[i_phase][i_data].mass_integral_target_global)){

                                balance_error_exists = 0;
                            }

                            i_while++;
                        }

                        if( (i_while == rcfd->Balance_Dict[i_phase][i_data].max_correction_loops)
                            && (fabs(  rcfd->Balance[i_phase][i_data].mass_integral_target_global
                                     - rcfd->Balance[i_phase][i_data].mass_integral_global)
                                < rcfd->Balance_Dict[i_phase][i_data].accuracy_level)){

                            report_warning0(rcfd, RCFD_REPORT_CATEGORY_RUN_BALANCE,
                                            "Balance for i_phase %d i_data %d beyond accuracy level %f",
                                            i_phase, i_data, rcfd->Balance_Dict[i_phase][i_data].accuracy_level);
                        }

                        break;
                    }
                    case RCFD_DATA_TEMPERATURE:
                    {
                        i_while = 0;

                        balance_error_exists = 1;

                        i_frame_c0 = rcfd->Rec->global_frame;
                        i_frame_c1 = rcfd->Rec->global_frame;

                        while((balance_error_exists) && (i_while < rcfd->Balance_Dict[i_phase][i_data].max_correction_loops)){

                            for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

                                _C(i_layer).data_swap[i_phase][i_cell][i_data] = 0.0;
                            }

                            available_exchange = 0.0;

                            for(i_face = 0; i_face < _FACE_DICT(i_layer).number_of_int_faces; i_face++){

                                c0 = _F(i_layer).c0[i_face];
                                c1 = _F(i_layer).c1[i_face];

                                if(_C(i_layer).data[i_phase][c0][i_data] > _C(i_layer).data[i_phase][c1][i_data]){

                                    c = c0; c0 = c1; c1 = c;
                                }

                                if(_C(i_layer).data[i_phase][c1][i_data] > _C(i_layer).data[i_phase][c0][i_data]){

#ifdef RCFD_LIB_ISLANDS
                                    i_frame_c0 = rcfd->Rec->global_frame_of_island[_C(i_layer).island_id[c0]];
                                    i_frame_c1 = rcfd->Rec->global_frame_of_island[_C(i_layer).island_id[c1]];
#endif

                                    if((_C(i_layer).volume[c0] * _C(i_layer).vof[i_frame_c0][c0][i_phase]) < (_C(i_layer).volume[c1] * _C(i_layer).vof[i_frame_c1][c1][i_phase])){

                                        vol_flip = _C(i_layer).volume[c0] * _C(i_layer).vof[i_frame_c0][c0][i_phase];
                                    }
                                    else{
                                        vol_flip = _C(i_layer).volume[c1] * _C(i_layer).vof[i_frame_c1][c1][i_phase];
                                    }

                                    available_exchange += 0.5
                                                        * (_C(i_layer).data[i_phase][c1][i_data] - _C(i_layer).data[i_phase][c0][i_data])
                                                        * vol_flip
                                                        * rcfd->Phase_Dict[i_phase].density
                                                        * rcfd->Phase_Dict[i_phase].heat_capacity
                                                        * rcfd->Solver_Dict->face_swap_max_per_loop;
                                }
                            }

                            if(rcfd->Balance_Dict[i_phase][i_data].type == RCFD_BALANCE_GLOBAL){

                                available_exchange = PRF_GRSUM1(available_exchange);

                                if(available_exchange > fabs(rcfd->Balance[i_phase][i_data].mass_error_global)){

                                    if(available_exchange > 0.0){

                                        exchange_ratio = fabs(rcfd->Balance[i_phase][i_data].mass_error_global)
                                                        / available_exchange;
                                    }
                                    else{

                                        exchange_ratio = 0.0;
                                    }
                                }
                                else{

                                    exchange_ratio = 1.0;
                                }
                            }
                            else{

                                exchange_ratio = 0.0;
                            }

                            for(i_face = 0; i_face < _FACE_DICT(i_layer).number_of_int_faces; i_face++){

                                c0 = _F(i_layer).c0[i_face];
                                c1 = _F(i_layer).c1[i_face];

                                if(_C(i_layer).data[i_phase][c0][i_data] > _C(i_layer).data[i_phase][c1][i_data]){

                                    c = c0, c1 = c0, c0 = c;
                                }

                                if(_C(i_layer).data[i_phase][c0][i_data] < _C(i_layer).data[i_phase][c1][i_data]){

#ifdef RCFD_LIB_ISLANDS
                                    i_frame_c0 = rcfd->Rec->global_frame_of_island[_C(i_layer).island_id[c0]];
                                    i_frame_c1 = rcfd->Rec->global_frame_of_island[_C(i_layer).island_id[c1]];
#endif

                                    if(   (_C(i_layer).volume[c0] * _C(i_layer).vof[i_frame_c0][c0][i_phase])
                                        < (_C(i_layer).volume[c1] * _C(i_layer).vof[i_frame_c1][c1][i_phase])){

                                        vol_flip = _C(i_layer).volume[c0] * _C(i_layer).vof[i_frame_c0][c0][i_phase];
                                    }
                                    else{

                                        vol_flip = _C(i_layer).volume[c1] * _C(i_layer).vof[i_frame_c1][c1][i_phase];
                                    }

                                    if(rcfd->Balance[i_phase][i_data].mass_error_global > 0.0){

                                        _C(i_layer).data_swap[i_phase][c0][i_data] += exchange_ratio
                                                                                    * (_C(i_layer).data[i_phase][c1][i_data] - _C(i_layer).data[i_phase][c0][i_data])
                                                                                    * vol_flip
                                                                                    * rcfd->Phase_Dict[i_phase].density
                                                                                    * rcfd->Phase_Dict[i_phase].heat_capacity
                                                                                    * rcfd->Solver_Dict->face_swap_max_per_loop;
                                    }

                                    if(rcfd->Balance[i_phase][i_data].mass_error_global < 0.0){

                                        _C(i_layer).data_swap[i_phase][c1][i_data] -= exchange_ratio
                                                                                    * (_C(i_layer).data[i_phase][c1][i_data] - _C(i_layer).data[i_phase][c0][i_data])
                                                                                    * vol_flip
                                                                                    * rcfd->Phase_Dict[i_phase].density
                                                                                    * rcfd->Phase_Dict[i_phase].heat_capacity
                                                                                    * rcfd->Solver_Dict->face_swap_max_per_loop;
                                    }
                                }
                            }

                            rcfd->Balance[i_phase][i_data].mass_integral = 0.0;

                            i_frame = rcfd->Rec->global_frame;

                            for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

#ifdef RCFD_LIB_ISLANDS
                                i_frame = rcfd->Rec->global_frame_of_island[_C(i_layer).island_id[i_cell]];
#endif

                                if((_C(i_layer).volume[i_cell]
                                        * _C(i_layer).vof[i_frame][i_cell][i_phase]
                                        * rcfd->Phase_Dict[i_phase].density
                                        * rcfd->Phase_Dict[i_phase].heat_capacity)
                                    > 0.0){

                                    _C(i_layer).data[i_phase][i_cell][i_data] = (
                                                                                    _C(i_layer).data[i_phase][i_cell][i_data]
                                                                                  * _C(i_layer).volume[i_cell]
                                                                                  * _C(i_layer).vof[i_frame][i_cell][i_phase]
                                                                                  * rcfd->Phase_Dict[i_phase].density
                                                                                  * rcfd->Phase_Dict[i_phase].heat_capacity
                                                                                  + _C(i_layer).data_swap[i_phase][i_cell][i_data]
                                                                                )
                                                                                /
                                                                                (
                                                                                    _C(i_layer).volume[i_cell]
                                                                                  * _C(i_layer).vof[i_frame][i_cell][i_phase]
                                                                                  * rcfd->Phase_Dict[i_phase].density
                                                                                  * rcfd->Phase_Dict[i_phase].heat_capacity
                                                                                );
                                }

                                rcfd->Balance[i_phase][i_data].mass_integral +=   _C(i_layer).data[i_phase][i_cell][i_data]
                                                                                * _C(i_layer).volume[i_cell]
                                                                                * _C(i_layer).vof[i_frame][i_cell][i_phase]
                                                                                * rcfd->Phase_Dict[i_phase].density
                                                                                * rcfd->Phase_Dict[i_phase].heat_capacity;
                            }

                            rcfd->Balance[i_phase][i_data].mass_integral_global = PRF_GRSUM1(rcfd->Balance[i_phase][i_data].mass_integral);
                            rcfd->Balance[i_phase][i_data].mass_error_global =    rcfd->Balance[i_phase][i_data].mass_integral_target_global
                                                                                - rcfd->Balance[i_phase][i_data].mass_integral_global;

                            if(fabs(rcfd->Balance[i_phase][i_data].mass_error_global) <
                                (rcfd->Balance_Dict[i_phase][i_data].accuracy_level * rcfd->Balance[i_phase][i_data].mass_integral_target_global)){

                                balance_error_exists = 0;
                            }

                            i_while++;
                        }

                        if(i_while == rcfd->Balance_Dict[i_phase][i_data].max_correction_loops){

                            report_warning0(rcfd, RCFD_REPORT_CATEGORY_RUN_BALANCE,
                                            "Balance for i_phase %d i_data %d beyond accuracy level %f",
                                            i_phase, i_data, rcfd->Balance_Dict[i_phase][i_data].accuracy_level);
                        }

                        break;
                    }
                }
            }
        }
#endif
    }

    void balance_control_conc_sum(rCFD_type* rcfd, const int i_phase)
    {
#if RP_NODE

#if 1   // local vars
        int     i_layer, i_data, i_cell;
        double  sum_of_conc;
#endif

        i_layer = rcfd->Solver->current_layer;

        if(rcfd->Solver_Dict->control_conc_sum_on){

            if(rcfd->Phase_Dict[i_phase].number_of_concentration_data > 1){

                for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

                    sum_of_conc = 0.0;

                    // crop conc values to [0..1], calc sum_of_conc
                    for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

                        if(rcfd->Data_Dict[i_phase][i_data].type == RCFD_DATA_CONCENTRATION){

                            if(_C(i_layer).data[i_phase][i_cell][i_data] < 0.0){

                                _C(i_layer).data[i_phase][i_cell][i_data] = 0.0;
                            }
                            else if(_C(i_layer).data[i_phase][i_cell][i_data] > 1.0){

                                _C(i_layer).data[i_phase][i_cell][i_data] = 1.0;
                            }

                            sum_of_conc += _C(i_layer).data[i_phase][i_cell][i_data];
                        }
                    }

                    // adapt conc values
                    for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

                        if(rcfd->Data_Dict[i_phase][i_data].type == RCFD_DATA_CONCENTRATION){

                            if(sum_of_conc > 0.0){

                                _C(i_layer).data[i_phase][i_cell][i_data] /= sum_of_conc;
                            }
                        }
                    }
                }
            }
        }
#endif
    }

#endif
