/*
 * Copyright (C) 2008-2012 Delft University of Technology - All rights reserved.
 *
 * Author: Reijer Idema
 */

/*
 * ForcingTerms class source.
 */
#include <iostream>			// USES cin
#include "ForcingTerms.h"

NONLINEARSOLVER_NAMESPACE_BEGIN

//============================= ABC ==========================================

#undef __FUNCT__
#define __FUNCT__ "ForcingTerms::ForcingTerms"
ForcingTerms::ForcingTerms(void) :
	mErr(0), mName(""), mRelativeTolerance(0), mAbsoluteTolerance(0),
	mDivergenceTolerance(0), mMaxIterationCount(0), mNonlinearAbsoluteTolerance(0)
{
} // ForcingTerms::ForcingTerms

#undef __FUNCT__
#define __FUNCT__ "ForcingTerms::~ForcingTerms"
/* virtual */ ForcingTerms::~ForcingTerms(void)
{ // Always provide a virtual destructor on base classes, even if it is empty.
} // ForcingTerms::~ForcingTerms

/*
#undef __FUNCT__
#define __FUNCT__ "ForcingTerms::SetLinearTolerances"
PetscErrorCode ForcingTerms::SetLinearTolerances(KSP ksp)
{
	mErr = KSPGetTolerances(ksp, &mRelativeTolerance, &mAbsoluteTolerance, &mDivergenceTolerance, &mMaxIterationCount); CHKERRQ(mErr);
	return 0;
} // ForcingTerms::SetLinearTolerances
*/

#undef __FUNCT__
#define __FUNCT__ "ForcingTerms::Initialize"
/* virtual */ PetscErrorCode ForcingTerms::Initialize(PetscReal relativeTolerance, PetscReal absoluteTolerance, PetscReal divergenceTolerance, PetscInt maxIterationCount, PetscReal nonlinearAbsoluteTolerance)
{
	mRelativeTolerance = relativeTolerance;
	mAbsoluteTolerance = absoluteTolerance;
	mDivergenceTolerance = divergenceTolerance;
	mMaxIterationCount = maxIterationCount;
	mNonlinearAbsoluteTolerance = nonlinearAbsoluteTolerance;
	return 0;
} // ForcingTerms::Initialize

/* virtual */ std::string ForcingTerms::GetName(void) const { return mName; }
/* virtual */ PetscReal ForcingTerms::GetRelativeTolerance(void) const { return mRelativeTolerance; }
/* virtual */ PetscReal ForcingTerms::GetAbsoluteTolerance(void) const { return mAbsoluteTolerance; }
/* virtual */ PetscReal ForcingTerms::GetDivergenceTolerance(void) const { return mDivergenceTolerance; }
/* virtual */ PetscInt ForcingTerms::GetMaxIterationCount(void) const { return mMaxIterationCount; }
/* virtual */ PetscReal ForcingTerms::GetNonlinearAbsoluteTolerance(void) const { return mNonlinearAbsoluteTolerance; }

//============================= FIXED ========================================

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsFixed::ForcingTermsFixed"
ForcingTermsFixed::ForcingTermsFixed(void) : mDefaultConvergedContext(PETSC_NULL)
{	// Base constructor is automatically called
	mName = "fixed";
} // ForcingTermsFixed::ForcingTermsFixed

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsFixed::~ForcingTermsFixed"
/* virtual */ ForcingTermsFixed::~ForcingTermsFixed(void)
{	// Base destructor is automatically called
	if (mDefaultConvergedContext != PETSC_NULL)
	{
		mErr = KSPDefaultConvergedDestroy(mDefaultConvergedContext); CHKERRV(mErr);
		mDefaultConvergedContext = PETSC_NULL;
	}
} // ForcingTermsFixed::~ForcingTermsFixed

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsFixed::Initialize"
/* virtual */ PetscErrorCode ForcingTermsFixed::Initialize(PetscReal relativeTolerance, PetscReal absoluteTolerance, PetscReal divergenceTolerance, PetscInt maxIterationCount, PetscReal nonlinearAbsoluteTolerance)
{
	if (mDefaultConvergedContext == PETSC_NULL)
	{
		mErr = KSPDefaultConvergedCreate(&mDefaultConvergedContext); CHKERRQ(mErr);
	}
	return ForcingTerms::Initialize(relativeTolerance, absoluteTolerance, divergenceTolerance, maxIterationCount, nonlinearAbsoluteTolerance);
} // ForcingTermsFixed::Initialize

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsFixed::TestConvergence"
/* virtual */ PetscErrorCode ForcingTermsFixed::TestConvergence(KSP ksp, PetscInt it, PetscReal rnorm, KSPConvergedReason* pReason)
{
	if (it == 0)
	{	// First linear iteration of a new Newton iteration, set tolerances
		mErr = PetscInfo4(PETSC_NULL, "Set tolerances: rtol=%e, atol=%e, dtol=%e, maxits=%d.\n", mRelativeTolerance, mAbsoluteTolerance, mDivergenceTolerance, mMaxIterationCount); CHKERRQ(mErr);
		mErr = KSPSetTolerances(ksp, mRelativeTolerance, mAbsoluteTolerance, mDivergenceTolerance, mMaxIterationCount); CHKERRQ(mErr);
	}
	// Check for convergence
	mErr = KSPDefaultConverged(ksp, it, rnorm, pReason, mDefaultConvergedContext); CHKERRQ(mErr);
	return 0;
} // ForcingTermsFixed::TestConvergence

//============================= DEMBO & STEIHAUG =============================

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsDemboSteihaug::ForcingTermsDemboSteihaug"
ForcingTermsDemboSteihaug::ForcingTermsDemboSteihaug(void)
	: mNonLinearIterationCount(-1), mNonLinearResidualNorm(0), mDefaultConvergedContext(PETSC_NULL)
{	// Base constructor is automatically called
	mName = "Dembo and Steihaug";
} // ForcingTermsDemboSteihaug::ForcingTermsDemboSteihaug

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsDemboSteihaug::~ForcingTermsDemboSteihaug"
/* virtual */ ForcingTermsDemboSteihaug::~ForcingTermsDemboSteihaug(void)
{	// Base destructor is automatically called
	if (mDefaultConvergedContext != PETSC_NULL)
	{
		mErr = KSPDefaultConvergedDestroy(mDefaultConvergedContext); CHKERRV(mErr);
		mDefaultConvergedContext = PETSC_NULL;
	}
} // ForcingTermsDemboSteihaug::~ForcingTermsDemboSteihaug

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsDemboSteihaug::Initialize"
/* virtual */ PetscErrorCode ForcingTermsDemboSteihaug::Initialize(PetscReal relativeTolerance, PetscReal absoluteTolerance, PetscReal divergenceTolerance, PetscInt maxIterationCount, PetscReal nonlinearAbsoluteTolerance)
{
	if (mDefaultConvergedContext == PETSC_NULL)
	{
		mErr = KSPDefaultConvergedCreate(&mDefaultConvergedContext); CHKERRQ(mErr);
	}
	mNonLinearIterationCount = -1;
	mNonLinearResidualNorm = 0;
	return ForcingTerms::Initialize(relativeTolerance, absoluteTolerance, divergenceTolerance, maxIterationCount, nonlinearAbsoluteTolerance);
} // ForcingTermsDemboSteihaug::Initialize

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsDemboSteihaug::TestConvergence"
/* virtual */ PetscErrorCode ForcingTermsDemboSteihaug::TestConvergence(KSP ksp, PetscInt it, PetscReal rnorm, KSPConvergedReason* pReason)
{
	if (it == 0)
	{
		// First linear iteration of a new Newton iteration
		++mNonLinearIterationCount;
		mNonLinearResidualNorm = rnorm;
		mErr = UpdateTolerances(); CHKERRQ(mErr);	// Calculate new tolerances
		mErr = PetscInfo4(PETSC_NULL, "Set tolerances: rtol=%e, atol=%e, dtol=%e, maxits=%d.\n", mRelativeTolerance, mAbsoluteTolerance, mDivergenceTolerance, mMaxIterationCount); CHKERRQ(mErr);
		mErr = KSPSetTolerances(ksp, mRelativeTolerance, mAbsoluteTolerance, mDivergenceTolerance, mMaxIterationCount); CHKERRQ(mErr);
	}
	// Check for convergence
	mErr = KSPDefaultConverged(ksp, it, rnorm, pReason, mDefaultConvergedContext); CHKERRQ(mErr);
	return 0;
} // ForcingTermsDemboSteihaug::TestConvergence

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsDemboSteihaug::UpdateTolerances"
/* virtual */ PetscErrorCode ForcingTermsDemboSteihaug::UpdateTolerances(void)
{
	// Relative tolerance is minimum of (1/2)^i and ||F_i||
	PetscReal power = 0.5;
	for (PetscInt i = 0; i < mNonLinearIterationCount; ++i) power /= 2;
	PetscReal rtol = PetscMin(power, mNonLinearResidualNorm);

	// Enforce lower bound for nonlinear residual norm to (comfortably) reach its absolute tolerance if there is no oversolving
	PetscReal rtolLowerBound = PetscMax(0.1*mNonlinearAbsoluteTolerance/mNonLinearResidualNorm, 1e-12);
	if (rtol < rtolLowerBound) rtol = rtolLowerBound;

	// Update relative tolerance
	mRelativeTolerance = rtol;
	return 0;
} // ForcingTermsDemboSteihaug::UpdateTolerances

//============================= EISENSTAT & WALKER ===========================

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsEisenstatWalker::ForcingTermsEisenstatWalker"
ForcingTermsEisenstatWalker::ForcingTermsEisenstatWalker(void) : mNonLinearIterationCount(-1),
	mLinearResidualNorm(0), mLinearResidualNormOld(0),
	mNonLinearResidualNorm(0), mNonLinearResidualNormOld(0),
	mRelativeToleranceOld(0), mDefaultConvergedContext(PETSC_NULL)
{	// Base constructor is automatically called
	mName = "Eisenstat and Walker";
} // ForcingTermsEisenstatWalker::ForcingTermsEisenstatWalker

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsEisenstatWalker::~ForcingTermsEisenstatWalker"
/* virtual */ ForcingTermsEisenstatWalker::~ForcingTermsEisenstatWalker(void)
{	// Base destructor is automatically called
	if (mDefaultConvergedContext != PETSC_NULL)
	{
		mErr = KSPDefaultConvergedDestroy(mDefaultConvergedContext); CHKERRV(mErr);
		mDefaultConvergedContext = PETSC_NULL;
	}
} // ForcingTermsEisenstatWalker::~ForcingTermsEisenstatWalker

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsEisenstatWalker::Initialize"
/* virtual */ PetscErrorCode ForcingTermsEisenstatWalker::Initialize(PetscReal relativeTolerance, PetscReal absoluteTolerance, PetscReal divergenceTolerance, PetscInt maxIterationCount, PetscReal nonlinearAbsoluteTolerance)
{
	if (mDefaultConvergedContext == PETSC_NULL)
	{
		mErr = KSPDefaultConvergedCreate(&mDefaultConvergedContext); CHKERRQ(mErr);
	}

	mNonLinearIterationCount = -1;
	mLinearResidualNorm = 0;
	mLinearResidualNormOld = 0;
	mNonLinearResidualNorm = 0;
	mNonLinearResidualNormOld = 0;
	mRelativeToleranceOld = 0;

	return ForcingTerms::Initialize(relativeTolerance, absoluteTolerance, divergenceTolerance, maxIterationCount, nonlinearAbsoluteTolerance);
} // ForcingTermsEisenstatWalker::Initialize

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsEisenstatWalker::TestConvergence"
/* virtual */ PetscErrorCode ForcingTermsEisenstatWalker::TestConvergence(KSP ksp, PetscInt it, PetscReal rnorm, KSPConvergedReason* pReason)
{
	if (it == 0)
	{
		// First linear iteration of a new Newton iteration, update Newton iteration data
		mLinearResidualNormOld = mLinearResidualNorm;
		mNonLinearResidualNormOld = mNonLinearResidualNorm;
		mNonLinearResidualNorm = rnorm;
		++mNonLinearIterationCount;
	}
	mLinearResidualNorm = rnorm;	// Save linear residual norm
	if (it == 0)
	{
		// First linear iteration of a new Newton iteration, set linear tolerances for this Newton iteration
		mRelativeToleranceOld = mRelativeTolerance;
		mErr = UpdateTolerances(); CHKERRQ(mErr);	// Calculate new tolerances
		mErr = PetscInfo4(PETSC_NULL, "Set tolerances: rtol=%e, atol=%e, dtol=%e, maxits=%d.\n", mRelativeTolerance, mAbsoluteTolerance, mDivergenceTolerance, mMaxIterationCount); CHKERRQ(mErr);
		mErr = KSPSetTolerances(ksp, mRelativeTolerance, mAbsoluteTolerance, mDivergenceTolerance, mMaxIterationCount); CHKERRQ(mErr);
	}

	// Check for convergence
	mErr = KSPDefaultConverged(ksp, it, rnorm, pReason, mDefaultConvergedContext); CHKERRQ(mErr);
	return 0;
} // ForcingTermsEisenstatWalker::TestConvergence

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsEisenstatWalker::UpdateTolerances"
/* virtual */ PetscErrorCode ForcingTermsEisenstatWalker::UpdateTolerances(void)
{
	PetscReal rtol = 0;
	if (mNonLinearIterationCount == 0)
	{
		rtol = 0.1;		// Relative tolerance in first nonlinear iteration
	}
	else if (mNonLinearResidualNorm >= mNonLinearResidualNormOld)
	{
		rtol = 0.1;		// Special tolerance in case of divergence, otherwise tolerance might get >= 1.
	}
	else
	{	// Calculate relative tolerance by the Eisenstat & Walker method
		rtol = fabs(mNonLinearResidualNorm - mLinearResidualNormOld) / mNonLinearResidualNormOld;
		PetscReal safeguard = pow(mRelativeToleranceOld, (1+sqrt(5))/2);
		if (safeguard > 0.1 && safeguard > rtol) rtol = safeguard;
	}

	// Enforce lower bound for nonlinear residual norm to (comfortably) reach its absolute tolerance if there is no oversolving
	PetscReal rtolLowerBound = PetscMax(0.1*mNonlinearAbsoluteTolerance/mNonLinearResidualNorm, 1e-12);
	if (rtol < rtolLowerBound) rtol = rtolLowerBound;

	// Update relative tolerance
	mRelativeTolerance = rtol;
	return 0;
} // ForcingTermsEisenstatWalker::UpdateTolerances

//============================= ITERATIONS ===================================

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsIterations::ForcingTermsIterations"
ForcingTermsIterations::ForcingTermsIterations(void) : mNonLinearIterationCount(-1)
{ // Base constructor is automatically called
	mName = "iterations";
} // ForcingTermsIterations::ForcingTermsIterations

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsIterations::~ForcingTermsIterations"
/* virtual */ ForcingTermsIterations::~ForcingTermsIterations(void)
{ // Base destructor is automatically called
} // ForcingTermsIterations::~ForcingTermsIterations

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsIterations::Initialize"
/* virtual */ PetscErrorCode ForcingTermsIterations::Initialize(PetscReal relativeTolerance, PetscReal absoluteTolerance, PetscReal divergenceTolerance, PetscInt maxIterationCount, PetscReal nonlinearAbsoluteTolerance)
{
	mNonLinearIterationCount = -1;
	return ForcingTerms::Initialize(relativeTolerance, absoluteTolerance, divergenceTolerance, maxIterationCount, nonlinearAbsoluteTolerance);
} // ForcingTermsIterations::Initialize

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsIterations::TestConvergence"
/* virtual */ PetscErrorCode ForcingTermsIterations::TestConvergence(KSP ksp, PetscInt it, PetscReal rnorm, KSPConvergedReason* pReason)
{
	if (it == 0)
	{
		// First linear iteration of a new Newton iteration
		++mNonLinearIterationCount;
		mErr = UpdateTolerances(); CHKERRQ(mErr);
		mErr = PetscInfo4(PETSC_NULL, "Set tolerances: rtol=%e, atol=%e, dtol=%e, maxits=%d.\n", mRelativeTolerance, mAbsoluteTolerance, mDivergenceTolerance, mMaxIterationCount); CHKERRQ(mErr);
		mErr = KSPSetTolerances(ksp, mRelativeTolerance, mAbsoluteTolerance, mDivergenceTolerance, mMaxIterationCount); CHKERRQ(mErr);
	}
	// Checking by iterations, do not check tolerances
	mErr = KSPSkipConverged(ksp, it, rnorm, pReason, NULL); CHKERRQ(mErr);
	return 0;
} // ForcingTermsIterations::TestConvergence

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsIterations::UpdateTolerances"
/* virtual */ PetscErrorCode ForcingTermsIterations::UpdateTolerances(void)
{
	// Ask user for desired iteration count...
	PetscErrorCode ierr = PetscPrintf(PETSC_COMM_WORLD, "Enter number of linear iterations for Newton iteration %i: ", mNonLinearIterationCount+1); CHKERRQ(ierr);
	PetscInt i;
	std::cin >> i;	// :TODO: FIND SOME OTHER FORM OF INPUT?
	mMaxIterationCount = i;
	return 0;
} // ForcingTermsIterations::UpdateTolerances

//============================================================================

NONLINEARSOLVER_NAMESPACE_END
