/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*ConvFail handle conergece failure */
#include "cddefines.h"
#include "cddrive.h"
#include "prt.h"
#include "phycon.h"
#include "pressure.h"
#include "dense.h"
#include "thermal.h"
#include "called.h"
#include "map.h"
#include "wind.h"
#include "conv.h"

/*ConvFail handle conergece failure */
void ConvFail(
	/* chMode is one of "pres", "eden", "ioni", "pops", "grai", "temp" */
	const char chMode[], /* chMode[5] */
	const char chDetail[] )
{
	double relerror ;

#	ifdef DEBUG_FUN
	fputs( "<+>ConvFail()\n", debug_fp );
#	endif

	/* >>chng 02 jun 15, add this branch */
	/* do not announce a convergence failure - this was an abort before
	 * convergence was acheived */
	if( lgAbort )
	{
		/* an abort is not one of the failures we deal with - simply return and
		 * let something else handle this */
		/*fprintf( ioQQQ, " ConvFail - abort has been set.\n");*/
#		ifdef DEBUG_FUN
		fputs( " <->ConvFail()\n", debug_fp );
#		endif
		return;
	}

	/* pressure failure */
	if( strcmp( chMode , "pres" )==0 )
	{
		/* record number of pressure failures */
		++conv.nPreFail;
		if( called.lgTalk )
		{
			fprintf( ioQQQ, 
				" PROBLEM  ConvFail %li, pressure not converged; itr %li, zone %.2f Te:%.3e Hden:%.4e curr Pres:%.4e Corr Pres:%.4e Pra/gas:%.3e\n", 
			  conv.nPreFail,
			  iteration,
			  fnzone, 
			  phycon.te, 
			  dense.gas_phase[ipHYDROGEN], 
			  pressure.PresTotlCurr, 
			  pressure.PresTotlCorrect,
			  pressure.pbeta);

			/* this identifies new dynamics that failed near the sonic point */
			if( fabs(pressure.PresGasCurr - pressure.PresRamCurr)/pressure.PresGasCurr < 0.1 &&
				((strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv < 0. ) )
			{
				fprintf( ioQQQ, 
					"\n PROBLEM continued, pressure not converged; we are stuck at the sonic point.\n\n");
				pressure.lgSonicPoint = TRUE;
			}
		}
	}

	/* electron density failure */
	else if( strcmp( chMode, "eden" ) == 0 )
	{
		/* record number of electron density failures */
		++conv.nNeFail;

		if( called.lgTalk )
		{
			fprintf( ioQQQ, 
				" PROBLEM  ConvFail %li, eden not converged itr %li zn %.2f correct=%.3e "
				"assumed=%.3e.", 
			  conv.nNeFail,
			  iteration ,
			  fnzone,
			  dense.EdenTrue, 
			  dense.eden
			  );

			/* some extra information that may be printed */
			/* heating cooling failure */
			if( !conv.lgConvTemp )
			{
				fprintf( ioQQQ, "  Temperature failure also." );
			}

			/* heating cooling failure */
			if( !conv.lgConvIoniz )
			{
				fprintf( ioQQQ, "  Ionization failure also." );
			}

			/* electron density is oscillating */
			if( conv.lgEdenOscl )
			{
				fprintf( ioQQQ, "  Electron density oscillating." );
			}

			/* heating cooling oscillating */
			if( conv.lgCmHOsc )
			{
				fprintf( ioQQQ, "  Heating-cooling oscillating." );
			}
		}
		fprintf( ioQQQ, " \n");
	}

	else if( strcmp( chMode, "ioni" ) == 0 )
	{
		/* ionization failure */
		++conv.nIonFail;
		if( called.lgTalk )
		{
			fprintf( ioQQQ, " PROBLEM  ConvFail %li, %s ionization not converged zone %.2f \n", 
			  conv.nIonFail, 
			  chDetail ,
			  fnzone  );
		}
	}

	else if( strcmp( chMode, "pops" ) == 0 )
	{
		/* populations failure */
		++conv.nPopFail;
		conv.lgConvPops = FALSE;
		if( called.lgTalk )
		{
			fprintf( ioQQQ, " PROBLEM  ConvFail %li, %s population not converged iteration %li zone=%.2f \n", 
			  conv.nPopFail, 
			  chDetail ,
			  iteration,
			  fnzone  );
		}
	}

	else if( strcmp( chMode, "grai" ) == 0 )
	{
		/* ionization failure */
		++conv.nGrainFail;
		if( called.lgTalk )
		{
			fprintf( ioQQQ, " PROBLEM  ConvFail %ld, a grain failure occurred iteration %li zone %.2f \n", 
			  conv.nGrainFail, 
			  iteration , 
			  fnzone );
		}
	}

	/* rest of routine is temperature failure */
	else if( strcmp( chMode, "temp" ) == 0 )
	{
		ASSERT( fabs((thermal.htot - thermal.ctot)/thermal.htot ) > conv.HeatCoolRelErrorAllowed );
		++conv.nTeFail;
		if( called.lgTalk )
		{
			fprintf( ioQQQ, 
				" PROBLEM  ConvFail %ld, Temp not converged itr %li zn %.2f Te=%.4e DTe=%.2e"
				" Htot=%.3e Ctot=%.3e rel err=%.3e rel tol:%.3e\n", 
			  conv.nTeFail, 
			  iteration ,
			  fnzone, 
			  phycon.te, 
			  thermal.dTemper ,
			  thermal.htot, 
			  thermal.ctot, 
			  (thermal.htot - thermal.ctot)/ thermal.htot,
			  conv.HeatCoolRelErrorAllowed );

			if( conv.lgTOscl && conv.lgCmHOsc )
			{
				fprintf( ioQQQ, " Temperature and d(Cool-Heat)/dT were both oscillating.\n" );
			}

			else if( conv.lgTOscl )
			{
				fprintf( ioQQQ, " Temperature was oscillating.\n" );
			}

			else if( conv.lgCmHOsc )
			{
				fprintf( ioQQQ, " d(Cool-Heat)/dT was oscillating.\n" );
			}

			/* not really a temperature failure, but something else */
			if( !conv.lgConvIoniz )
			{
				fprintf( ioQQQ, " Solution not converged due to %10.10s\n", 
				  conv.chConvIoniz );
			}
		}
	}
	else
	{
		fprintf( ioQQQ, " ConvFail called with insane mode %s detail %s\n", 
		  chMode , chDetail );
		ShowMe();
		puts( "[Stop in ConvFail]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* increment total number of failures */
	++conv.nTotalFailures;

	/* now see how many total failures we have, and if it is time to abort */
	/* remember which zone this is */
	conv.ifailz[MIN2(conv.nTotalFailures,10)-1] = nzone;

	/* remember the relative error
	 * convert to single precision for following max, abs (vax failed here) */
	relerror = fabs((thermal.htot-thermal.ctot)/ thermal.htot);

	conv.failmx = MAX2(conv.failmx,(float)relerror);
	/*conv.failmx = (float)fmax(conv.failmx,fabs((float)((thermal.htot-thermal.ctot)/
	  thermal.htot)));*/

	/* this branch is non-abort exit - we have not exceeded the limit to the number of failures */
	if( conv.nTotalFailures < conv.LimFail )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->ConvFail()\n", debug_fp );
#		endif
		return;
	}

	fprintf( ioQQQ, " Stop due to excessive convergence failures - there have been %ld so far. \n", 
		conv.LimFail );
	fprintf( ioQQQ, " This limit can be reset with the FAILURES command.\n" );
	fprintf( ioQQQ, "\n The input commands that lead to this failure are the following:\n" );
	/* now output the commands */
	cdPrintCommands(ioQQQ);
	fprintf( ioQQQ, "\n" );

	/* if due to pressure failures then recommend looking at pressure map */
	if( conv.nPreFail==conv.nTotalFailures )
	{
		fprintf( ioQQQ, " These were all pressure failures - we may be near an unstable point in the cooling curve. \n");
		fprintf( ioQQQ, " The PUNCH PRESSURE HISTORY command will show the n-T-P curve, and may be interesting.\n\n");
	}

	/* punt */
	if( conv.lgMap )
	{
		/* only do map if requested */
		/* adjust range of punting map */
		map.RangeMap[0] = (float)(phycon.te/100.);
		map.RangeMap[1] = (float)MIN2(phycon.te*100.,9e9);
		/* need to make printout out now, before disturbing soln with map */
		PrtZone();
		map.lgMapBeingDone = TRUE;
		map_do(ioQQQ,"punt");
	}

	/* return out from here and let iter_end_check catch lgAbort set,
	 * and generate normal output there */
	lgAbort = TRUE;
	if( called.lgTalk )
	{
		fprintf( ioQQQ, " ConvFail sets lgAbort since nTotalFailures=%ld is >= LimFail=%ld\n", 
		  conv.nTotalFailures, 
		  conv.LimFail );
		fprintf( ioQQQ, " This limit can be reset with the FAILURES command.\n");
	}

#	ifdef DEBUG_FUN
	fputs( " <->ConvFail()\n", debug_fp );
#	endif

	return;
}
