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

#ifndef RCFD_RUN
#define RCFD_RUN

#include "rCFD_layer.h"
#include "rCFD_post.h"
#include "rCFD_balance.h"


    void run_update_state(rCFD_type* rcfd)
    {
#if RP_NODE
        rcfd->Solver->current_state = 0;

        rcfd->Solver->upcoming_state = 0;
#endif // RP_NODE
    }

    void run_update_layer(rCFD_type* rcfd)
    {
        int prev_layer = rcfd->Solver->current_layer;

        rcfd->Solver->current_layer = rcfd->rCFD_UDF->_user_set_layer(rcfd, rcfd->Solver->current_layer);

        if(rcfd->Solver->current_layer != prev_layer){

            layer_map_from_to(rcfd, prev_layer, rcfd->Solver->current_layer);
        }
    }

/**
 * @brief helper fct for run_update_frame: normalizes (user set) frame_id
 * @param rcfd pointer to previously created rCFD instance
 */
 void normalize_frame_id(rCFD_type* rcfd){
#if RP_NODE

	if(myid == 0){
		
		int    i_id;

		int    frame_id_length = rcfd->Rec_Dict->frame_id_vector_length;
		
		double *frame_id =       rcfd->Rec->frame_id;
		double *frame_id_means = rcfd->Rec->frame_id_means; 
		double *frame_id_vars =  rcfd->Rec->frame_id_vars;
		
		for(i_id = 0; i_id < frame_id_length; i_id++){
		
			frame_id[i_id] -= frame_id_means[i_id];
			
			if(frame_id_vars[i_id] > 0.0){
				
				frame_id[i_id] /= frame_id_vars[i_id];
			}
		}

		report_info0(rcfd, 
					 RCFD_REPORT_CATEGORY_READ_DB, 
					 "run_update_frame: normalized frame_id");					
	}
#endif	
}

/**
 * @brief helper fct for run_update_frame: gets new frame id by comparing frame_id to DB_frame_id
 *     major code options:    frame_id_avoid_immediate_seq_repetition_on
 * @param rcfd pointer to previously created rCFD instance
 * @param frame_id current frame identification vector
 * @return index of starting frame of new sequence or 0 if error
 */
int get_global_frame_by_frame_id(rCFD_type* rcfd, double *frame_id)
{
#if RP_NODE

    if(myid != 0) return 0;

    int    i_frame, i_id;

    int    frame_count =         rcfd->Solver_Dict->number_of_frames;
    int    frame_id_length =     rcfd->Rec_Dict->frame_id_vector_length;
    int    frame_id_seq_length = rcfd->Rec_Dict->frame_id_sequence_length;

    int    min_dist_id =  0;
    int    new_frame_id = 0;
    double dist =        0.0;
    double min_dist =    99999.9;

    double **DB_frame_id = rcfd->Rec->frame_id_database;

    // get min_dist_id
    for(i_frame = 0; i_frame < (frame_count - frame_id_seq_length); i_frame++){

        dist = 0.0;

        for(i_id = 0; i_id < frame_id_length; i_id++){

            dist += fabs(frame_id[i_id] - DB_frame_id[i_frame][i_id]);
        }

        if(min_dist > dist){

            min_dist =    dist;
            min_dist_id = i_frame;
        }
    }

    new_frame_id = min_dist_id - (min_dist_id % frame_id_seq_length);
    
    // avoid immediate replay
    #define frame_id_avoid_immediate_seq_repetition_on  1   // should be var in Rec_Dict
    if(frame_id_avoid_immediate_seq_repetition_on){
        
        if(new_frame_id == (rcfd->Rec->global_frame - frame_id_seq_length + 1)){    
            
            new_frame_id += frame_id_seq_length;
            
            if(new_frame_id > frame_count){
            
                new_frame_id -= 2 * frame_id_seq_length;
            }
            
            if(new_frame_id < 0){
            
                new_frame_id = 0;
            }
        }
    }
    
    report_debug0(rcfd, RCFD_REPORT_CATEGORY_RUN_UPDATE,
           "run_update_frame: get_global_frame_by_frame_id: new frame %d, min_dist %e",
           new_frame_id, min_dist);
           
    return new_frame_id;

#else
    return 0;
#endif
}
   
/**
 * @brief sets the next frame of c2c shifts
 *     major code options:    recurrence_process_on, frame_id_on, adapt_vof_stitching_on
 *     existing limitations:  timestep width has to be adapted for multi-layer
 * @param rcfd pointer to previously created rCFD instance
 */
void run_update_frame(rCFD_type* rcfd)
{

// TODO at the moment run_update_frame works only for single layer

#if 1   // local vars
    int     i_state, i_state2, i_layer, i_rec;
    int     i_rec_max;
#ifdef RCFD_LIB_ISLANDS
    int     i_island
#endif
#if RP_HOST
    double  rand_real;
#endif // RP_HOST
#endif

    i_state     = rcfd->Solver->current_state;
    i_state2    = rcfd->Solver->upcoming_state;
    i_layer     = rcfd->Solver->current_layer;

    rcfd->Rec->jumped_at_last_frame = 0;

#ifndef RCFD_LIB_ISLANDS

    rcfd->Rec->prev_global_frame = rcfd->Rec->global_frame;

    if (rcfd->Rec->global_frame == (rcfd->Solver_Dict->number_of_frames - 1)){

        rcfd->Rec->jumped_at_last_frame = 1;
    }
#else
    for(i_island = 0; i_island < rcfd->Solver_Dict->number_of_islands; i_island++){

        rcfd->Rec->prev_global_frame_of_island[i_island] = rcfd->Rec->global_frame_of_island[i_island];

        if (rcfd->Rec->global_frame_of_island[i_island] == (rcfd->Solver_Dict->number_of_frames - 1)){

            rcfd->Rec->jumped_at_last_frame = 1;
        }
    }
#endif

    if(rcfd->Solver_Dict->recurrence_process_on){
        
        report_debug0(rcfd, RCFD_REPORT_CATEGORY_RUN_UPDATE,
               "[%d] run_update_frame: start, i_frame %d, rcfd->Rec->frame_in_sequence %d, rcfd->Rec->sequence_length %d",
               myid, rcfd->Rec->global_frame, rcfd->Rec->frame_in_sequence, rcfd->Rec->sequence_length);

#if RP_HOST // init rand. generator
        
#ifdef RCFD_STANDALONE
        if(myid == 0)
#endif
        {
            if(rcfd->Solver->global_run_counter == 0){

                unsigned int seed;

                if(rcfd->Rec_Dict->seed < 0){

                    seed = (unsigned int)time(0);
                }
                else{

                    seed = (unsigned int)rcfd->Rec_Dict->seed;
                }

                srand(seed);

#ifdef RCFD_STANDALONE
                report_print0(rcfd, "   --> Initialized random generator with seed %u", seed);
#else
                Message("\n\nInitialized random generator with seed %u", seed);
#endif
            }
        }
#endif // RP_HOST

        i_rec_max = 1;

        // TODO While the general multilayer strategy seems to be ok, it has to be refined for
        //      timestep widths unequal to Dt_1 = 2* Dt_0 ...
        for(i_rec = 0; i_rec < i_layer; i_rec++){

            i_rec_max *= 2;
        }
        for(i_rec = 0; i_rec < i_rec_max; i_rec++){

            if(rcfd->Rec->frame_in_sequence < (rcfd->Rec->sequence_length - 1)){

                rcfd->Rec->frame_in_sequence++;
                
                if(rcfd->Rec_Dict->format != RCFD_REC_FRAME_ID){

                    /* set next frame (within sequence) */
#ifndef RCFD_LIB_ISLANDS

                    if(rcfd->Rec->global_frame < (rcfd->Solver_Dict->number_of_frames - 1)){

                        rcfd->Rec->global_frame++;
                    }
                    else{

                        rcfd->Rec->global_frame = rcfd->Rec->jumps[i_state][i_state2][0][(rcfd->Solver_Dict->number_of_frames-1)];
                    }
#else
                    for(i_island = 0; i_island < rcfd->Solver_Dict->number_of_islands; i_island++){

                        if(rcfd->Rec->global_frame_of_island[i_island] < (rcfd->Solver_Dict->number_of_frames - 1)){

                            rcfd->Rec->global_frame_of_island[i_island]++;
                        }
                        else{

                            rcfd->Rec->global_frame_of_island[i_island] = rcfd->Rec->jumps[i_state][i_state2][i_island][(rcfd->Solver_Dict->number_of_frames-1)];
                        }
                    }
#endif
                }
                else if (rcfd->Rec_Dict->format == RCFD_REC_FRAME_ID){
                    
                    rcfd->Rec->global_frame++; // cannot exceed frame_count
                }
            }
            else{ /* new rCFD_seq */

                rcfd->Rec->frame_in_sequence = 0;
                
                if(rcfd->Rec_Dict->format != RCFD_REC_FRAME_ID){
#if RP_HOST
#ifdef RCFD_STANDALONE
                    if(myid == 0)
#endif // RCFD_STANDALONE
                    {
                        rand_real = (double)rand()/(double)RAND_MAX;

                        rcfd->Rec->sequence_length = rcfd->Rec_Dict->min_seq_length + (int)(rand_real*(double)(rcfd->Rec_Dict->max_seq_length - rcfd->Rec_Dict->min_seq_length));
                    }
#endif // RP_HOST
                    host_to_node_int_1(rcfd->Rec->sequence_length);


                    /* set next frame (new sequence) */
#ifndef RCFD_LIB_ISLANDS

                    rcfd->Rec->global_frame = rcfd->Rec->jumps[i_state][i_state2][0][rcfd->Rec->global_frame];
                    
                    report_debug0(rcfd, RCFD_REPORT_CATEGORY_RUN_UPDATE,
                           "[%d] run_update_frame: rcfd->Rec_Dict->format %d, new frame %d",
                           myid, rcfd->Rec_Dict->format, rcfd->Rec->global_frame);

#else
                    for(i_island = 0; i_island < rcfd->Solver_Dict->number_of_islands; i_island++){

                        rcfd->Rec->global_frame_of_island[i_island] = rcfd->Rec->jumps[i_state][i_state2][i_island][rcfd->Rec->global_frame_of_island[i_island]];
                    }
#endif
                }
                else if (rcfd->Rec_Dict->format == RCFD_REC_FRAME_ID){
#if RP_NODE
                    double *frame_id = rcfd->Rec->frame_id;

                    rcfd->rCFD_UDF->_user_set_frame_id(rcfd, frame_id);

					if( rcfd->Rec_Dict->frame_id_normalize_on){
						
						normalize_frame_id(rcfd);
					}

                    if(myid == 0){

                        rcfd->Rec->global_frame = get_global_frame_by_frame_id(rcfd, frame_id);
                    }
#endif 
                    PRF_BCAST_INT(&rcfd->Rec->global_frame, 1, 0); 
                }
            }
        }
    }
    else{

#ifndef RCFD_LIB_ISLANDS

        rcfd->Rec->global_frame = 0;
#else
        for(i_island = 0; i_island < rcfd->Solver_Dict->number_of_islands; i_island++){

            rcfd->Rec->global_frame_of_island[i_island] = 0;
        }
#endif
    }

    /* eventually, redistribute data on stitched vof fields */
    if(rcfd->Rec_Dict->adapt_vof_stitching_on){
#if RP_NODE
        layer_map_data_to_new_vof(rcfd);
#endif // RP_NODE
    }

	PRF_GSYNC();
}

    void run_update_balance(rCFD_type* rcfd)
    {
#if RP_NODE
        int     i_layer, i_phase;

        i_layer = rcfd->Solver->current_layer;

        for(i_phase = 0; i_phase < rcfd->Solver_Dict->number_of_phases; i_phase++){

            balance_init_source(rcfd, i_layer, i_phase);
        }
#endif // RP_NODE
    }

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

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

#endif

        i_layer = rcfd->Solver->current_layer;

        // init _C(i_layer).weight_after_shift/swap and _C(i_layer).data_shift/swap
        for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

            _C(i_layer).weight_after_shift[i_cell] =  0.0;
            _C(i_layer).weight_after_swap[i_cell] =   0.0;

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

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

        // init balance terms
        {
            balance_init_cross_partition_fluxes(rcfd, i_layer, i_phase);
        }

#endif // RP_NODE
    }

#endif
