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

/*
 * ForcingTerms class source.
 */
#include "ForcingTerms.h"


#include <iostream>


POWERFLOWER_NAMESPACE_BEGIN

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

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

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsTest::Initialize"
/* virtual */ PetscErrorCode ForcingTermsTest::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);
} // ForcingTermsTest::Initialize

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsTest::TestConvergence"
/* virtual */ PetscErrorCode ForcingTermsTest::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;
} // ForcingTermsTest::TestConvergence

#undef __FUNCT__
#define __FUNCT__ "ForcingTermsTest::UpdateTolerances"
/* virtual */ PetscErrorCode ForcingTermsTest::UpdateTolerances(void)
{
	PetscReal baseInitialResidualNorm = 5122.435401465;
	PetscReal x, y, xThreshold = 3, yThreshold = -2.21;

	// Get current x
	y = log10(mNonLinearResidualNorm/baseInitialResidualNorm);
	if (y >= yThreshold)
	{
		// y = a*x, a = yThreshold/xThreshold
		x = (xThreshold/yThreshold)*y;
	}
	else
	{
		// y = yThreshold + a*(x-3)^2, a = -0.5976088777
		x = 3 + sqrt((y-yThreshold)/-0.5976088777);
	}

	// Estimate norm for next x
	x = x + 1;
	if (x <= xThreshold)
	{
		// y = a*x, a = yThreshold/xThreshold
		y = (yThreshold/xThreshold)*x;
	}
	else
	{
		// y = yThreshold + a*(x-3)^2, a = -0.5976088777
		y = yThreshold - 0.5976088777 * (x-3) * (x-3);
	}
	PetscReal estimatedResidualNorm = baseInitialResidualNorm*pow(10, y);
	mRelativeTolerance = estimatedResidualNorm / mNonLinearResidualNorm;

/*
	if (mNonLinearIterationCount == 0)
		mRelativeTolerance = 0.02;
	else if (mNonLinearIterationCount == 1)
		mRelativeTolerance = 0.003;
	else if (mNonLinearIterationCount == 2)
		mRelativeTolerance = 0.0004;
	else if (mNonLinearIterationCount == 3)
		mRelativeTolerance = 0.00005;
	else
		mRelativeTolerance = 0.00001;
*/

	// Enforce lower bound for nonlinear residual norm to (comfortably) reach its absolute tolerance if there is no oversolving
	PetscReal rtolLowerBound = 0.1*mNonlinearAbsoluteTolerance/mNonLinearResidualNorm;
	if (mRelativeTolerance < rtolLowerBound) mRelativeTolerance = rtolLowerBound;

	return 0;
} // ForcingTermsTest::UpdateTolerances

POWERFLOWER_NAMESPACE_END
