/**
* @file    main.cpp
* @author  Sanbai Li
* @date
* @brief   KarstVolSys simulator framework
* @version
* Copyright (c) 2019, NEU & SINOPEC.
*/

#include "DataInput.h"
#include "DataOutput.h"
#include "SimHM.h"

#include "DDM_2D.h"

using namespace std;
using namespace KarstVolSys;

int main(int argc, char* argv[]) {
	using namespace KarstVolSys;
	PrintVersion("2019 ", cout);
	char filename[MAX_STRING_LEN];
	bool if_echo = true;

	if (argc == 1) TypeInName(filename);
	else strcpy_s(filename, argv[1] + 1);

	if (argc == 3) if_echo = (strcmp(argv[2] + 1, "echo_on") == 0);

	DataBase AllData;
	FileReader filereader(if_echo);
	
	bool foundfile = filereader.InputFile(filename, &AllData); 
	
	if (foundfile){

		AllData.DATAPROCESS();
		
		SimHM MODEL(&AllData);

		MODEL.Simulate();

	}
	if (if_echo){
		cout << "Press \"Enter\" to quit." << endl;
		if (argc == 1)getchar();
		getchar();
	}
	return 0;
}

void SimHM::Simulate(){
	//
	//TerzaghiMdl();
	//MandelMdl();
	//THMMdl();

	int	start = clock(), t_end;
	int TmStepCounter = 0;
	int RealTSCount = 0;
	int NtCount = 1, Total_NtCount = 0;
	int CutBackCount = 0;

	int maxTmStps = AllData->maxTimeSteps;
	int maxIteNum = AllData->max_itr;
	double TimeSim = AllData->T_start;
	double End_time = AllData->T_end;
	double MindT = AllData->min_dt;
	double MaxdT = AllData->max_dt;
	double dT = AllData->dT;
	double Current_dT = dT;

	bool isConvg = false;
	bool dx_acpt = true;
	bool X_acpt = true;
	bool isRollback = false;
	bool isFractured = false;
	bool isNewlyFractured = false;
	
	Init(AllData->isRestart); //Initilize flow and geomechanical parts; Allow for restart option
	
	Open_OUTPUT_DATA_VaryVars(); //Open files to write

	while (true) {
 		int t_start = clock();
		// **** Newton iterations ****
		for (NtCount = 1;; ++NtCount, ++Total_NtCount){

			// ****The error is small enough****
			if (dx_acpt && X_acpt && NtCount > 1){
				--NtCount;
				break;
			}

			// ****If Newton iteration fails****
			if (NtCount > maxIteNum) {
				break;
			}

			// ****Compute properties and flows, build matrix****
			Active(NtCount, TimeSim, dT, isFractured); //Modify the elastic matrix and perm. in function ActiveSolid()

			// **** If not fail, solve the linear system
			SolveEQs(Total_NtCount, isFractured, CSR_NNZ_); //isFractured = true means that reconstructing the Jacbian Matrix is needed!

			// Check the primary variables are fall in the tolerant range
            dx_acpt = CheckDx(dT);

			// Update primary and sencondary variables for next iteration
			X_acpt = UpdateX(dT);
			UpdatePerm4EDFM();   //Karst evolution, added in 2018/7
		}
		TimeSim += dT;
		TmStepCounter++;

		if (AllData->isFracing){
			if (AllData->hasNFs)isNewlyFractured = FailureCheckinNFR();  // Hydraulic fracturing in naturally fractured reservoirs
			else isNewlyFractured = FailureCheck();                      // Hydraulic fracturing in intact reservoirs
			if (AllData->isEDFM) isNewlyFractured = FailureCheck4EDFM(); /* 1. Stress state checking for emdedded natural fractures
																		    2. Update fracture permeability for emdedded fracture segments */
		} else;

		if (!isFractured) isFractured = isNewlyFractured;
		bool isMinDt = abs(dT - MindT)<1.0e-2 || dT <= MindT;
		isRollback = ((isNewlyFractured || NtCount > maxIteNum) && !isMinDt);

		//if (false/*isRollback*/){
		if (isRollback){
			TimeSim -= dT;
			RollBackX0();
			++CutBackCount;
			cout << "CutBackCount = " << CutBackCount;
		}
		else {
			cout << "/**************************************/" << endl;
			cout << "  Reach convergence at N=" << TmStepCounter << " time step!\n";
			cout << "/**************************************/" << endl; //Sleep(2 * 100);
			UpdateX0(dT);

			//
			//UpdatePerm4EDFM(); //Karst evolution, added in 2018/7
			AllData->UpdateTransmisibility();
			//
			
			/*
			  Output for 2D or 3D visualization
			*/
			Write_OUTPUT_DATA_VaryVars(TimeSim, dT); //Include data needed for restarting
			OUTPUT_DATA_StateVars(TimeSim);
		}
		/**************************************************************************/
		if (TimeSim >= End_time || TmStepCounter >= maxTmStps) {
			break;
		}
		//Get the size of timestep 
		SET_NEXT_STEP(TimeSim, dT, MaxdT, MindT, isRollback, Current_dT, NtCount);

		//Elapsed CPU time
		total_time = (double)(clock() - start);

		if (isRollback){
			Wasted_time += (double)(clock() - t_start);
		}
	}

	if (!AllData->isFracing) OUTPUT_FracForm << "Non-fracturing case. No info available!" << endl;

	ostream::fmtflags log_old = cout.flags();
	KarstVolSysLog << setiosflags(ios::fixed) << setprecision(2);
	KarstVolSysLog << "-------------------- Performance Statistics ---------------------------\n";
	KarstVolSysLog << "***********************************************************\n";
	KarstVolSysLog << "Performed:  " << TmStepCounter << " time steps (" << CutBackCount << " wasted)" << "\n";
	KarstVolSysLog << "            " << Total_NtCount << " Newton steps \n";
	KarstVolSysLog << "***********************************************************\n";
	KarstVolSysLog << "Total modeling time:                      " << setw(14) << TimeSim << " s\n";
	//KarstVolSysLog << "The current timestep:                     " << setw(14) << Current_dT << " s\n";
	KarstVolSysLog << "Total CPU time:                           " << setw(14) << total_time << " ms\n";
	KarstVolSysLog << "Wasted CPU time:                          " << setw(14) << Wasted_time << " ms\n";
	KarstVolSysLog << "Initialization time:                      " << setw(14) << InitTime << " ms\n";
	KarstVolSysLog << "Solver time:                              " << setw(14) << solver_time << " ms\n";
	KarstVolSysLog << "UpdatePerm2ReconJac_time:                 " << setw(14) << UpdatePerm2ReconJac_time << " ms\n";
	KarstVolSysLog << "ActiveSolid:                              " << setw(14) << ActiveSolid_time << " ms\n";
	KarstVolSysLog << "ActiveFlow:                               " << setw(14) << ActiveFlow_time << " ms\n";
	KarstVolSysLog << "ActiveOthers:                             " << setw(14) << ActiveOthers_time << " ms\n";
	KarstVolSysLog << "Check Convengence or failure state time:  " << setw(14) << CheckDvar2isFailure_time << " ms\n";
	KarstVolSysLog << "WriteData time:                           " << setw(14) << WriteData_time << " ms\n";
	KarstVolSysLog << "***********************************************************\n";
	KarstVolSysLog << "\----------------  Simulation of case \'" << AllData->casename << "\' complete -------------------\n" << endl;

	KarstVolSysLog.flags(log_old);
	KarstVolSysLog << setprecision(6);
	Close_OUTPUT_DATA_VaryVars();
}