/*!\file: transient_3d_core.cpp
 * \brief: core of the transient_3d solution 
 */ 

#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "../toolkits/toolkits.h"
#include "../objects/objects.h"
#include "../shared/shared.h"
#include "../io/io.h"
#include "../EnumDefinitions/EnumDefinitions.h"
#include "./solutions.h"
#include "../modules/modules.h"
#include "../include/include.h"
#include "../solvers/solvers.h"
#include <float.h>

void transient_core(FemModel* femmodel){

	/*parameters: */
	double finaltime,dt,yts;
	bool   control_analysis,isdiagnostic,isprognostic,isthermal,isgroundingline,isenthalpy;
	bool   dakota_analysis=false;
	bool   time_adapt=false;
	int    solution_type;
	int    output_frequency;
	int    dim,groundingline_migration;
	int    numoutputs         = 0;
	int    *requested_outputs = NULL;
	
	/*intermediary: */
	int    step;
	double time;

	//first recover parameters common to all solutions
	femmodel->parameters->FindParam(&dim,MeshDimensionEnum);
	femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
	femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
	femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
	femmodel->parameters->FindParam(&control_analysis,InversionIscontrolEnum);
	femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
	femmodel->parameters->FindParam(&output_frequency,SettingsOutputFrequencyEnum);
	femmodel->parameters->FindParam(&time_adapt,TimesteppingTimeAdaptEnum);
	femmodel->parameters->FindParam(&isdiagnostic,TransientIsdiagnosticEnum);
	femmodel->parameters->FindParam(&isprognostic,TransientIsprognosticEnum);
	femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum);
	femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
	femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
	if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
	femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
	if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,TransientRequestedOutputsEnum);

	/*initialize: */
	step=0;
	time=0;

	/*for qmu analysis, reinitialize velocity so that fake sensitivities do not show up as a result of a different restart of the convergence at each trial.*/
	if(dakota_analysis){
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVxEnum,VxEnum);
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVyEnum,VyEnum);
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVzEnum,VzEnum);
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuPressureEnum,PressureEnum);
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuThicknessEnum,ThicknessEnum);
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuSurfaceEnum,SurfaceEnum);
		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuBedEnum,BedEnum);
		if(isthermal && dim==3){
			InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuTemperatureEnum,TemperatureEnum);
			InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuMeltingEnum,BasalforcingsMeltingRateEnum);
		}
	}

	while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime.

		/*Increment*/
		if(time_adapt){
			TimeAdaptx(&dt,femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); 
			if(time+dt>finaltime) dt=finaltime-time;
			femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);
		}
		step+=1;
		time+=dt;
		femmodel->parameters->SetParam(time,TimeEnum);

		_printf_(VerboseSolution(),"iteration %i/%g  time [yr]: %-7.3g (time step: %.2g)\n",step,floor((finaltime-time)/dt)+step,time/yts,dt/yts);

		if(isthermal && dim==3){
			_printf_(VerboseSolution(),"   computing temperatures:\n");
			#ifdef _HAVE_THERMAL_
			if(isenthalpy==0){
				thermal_core_step(femmodel,step,time);
			}
			else{
				enthalpy_core_step(femmodel,step,time);
			}
			#else
			_error_("ISSM was not compiled with thermal capabilities. Exiting");
			#endif
		}

		if(isdiagnostic){
			_printf_(VerboseSolution(),"   computing new velocity\n");
			#ifdef _HAVE_DIAGNOSTIC_
			diagnostic_core(femmodel);
			#else
			_error_("ISSM was not compiled with diagnostic capabilities. Exiting");
			#endif
		}

		if(isprognostic){
			_printf_(VerboseSolution(),"   computing new thickness\n");
			prognostic_core(femmodel);
			_printf_(VerboseSolution(),"   updating vertices positions\n");
			UpdateVertexPositionsx(femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); 
		}

		if(isgroundingline){
			_printf_(VerboseSolution(),"   computing new grounding line position\n");
			#ifdef _HAVE_GROUNDINGLINE_
			GroundinglineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
			#else
			_error_("ISSM was not compiled with grounding line migration capabilities. Exiting");
			#endif
		}

		//save should be done by the analysis... see control_core
		if(solution_type==TransientSolutionEnum && !control_analysis && (step%output_frequency==0 || time==finaltime)){
			_printf_(VerboseSolution(),"   saving results\n");
			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzMeshEnum,step,time);
			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time);
			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time);
			if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,step,time);
			if(isdiagnostic)InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VelEnum,step,time);
			if(isdiagnostic)InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time);
			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time);
			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time);
			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time);
			if(dim==3 && isthermal) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,step,time);
			if(isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,WaterfractionEnum,step,time);
			if(isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EnthalpyEnum,step,time);
			if(!isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalforcingsMeltingRateEnum,step,time);
			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceforcingsMassBalanceEnum,step,time);
			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaskElementonfloatingiceEnum,step,time);
			RequestedOutputsx(femmodel->results,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,requested_outputs,numoutputs,step,time);

			/*unload results*/
			_printf_(VerboseSolution(),"   saving temporary results\n");
			OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,&femmodel->results);
		}
	}

	/*Free ressources:*/
	xfree((void**)&requested_outputs);
}
