#include "Sundance.h"

/** \example inlineRadDiffusion1D.cpp

Solve the radiation diffusion equation laplacian(T^4)=0 on the unit line.
This is a nonlinear problem, and is solved with a backtracking Newton's method.
In this example, the BT Newton algorithm is hardwired in main().

The exact solution is T(x)=[T(0)^4 + (T(1)^4-T(0)^4) x]^(1/4).

 */


int main(int argc, void** argv)
{
  try
    {
      Sundance::init(&argc, &argv);

      /* create a simple mesh on the unit line */
      int n = 40;
      MeshGenerator mesher = new LineMesher(0.0, 1.0, n);
      Mesh mesh = mesher.getMesh();

      /* define an expression representing the x-coordinate function */
      Expr x = new CoordExpr(0);

      /* define boundary cell sets */
      CellSet boundary = new BoundaryCellSet();
      CellSet left = boundary.subset( fabs(x - 0.0) < 1.0e-10 );
      CellSet right = boundary.subset( fabs(x - 1.0) < 1.0e-10 );

      /* varT is the variation of temperature */
      Expr varT = new TestFunction(new Lagrange(2));
      /* dT is the newton step; this is the unknown in the linearized problem. */
      Expr dT = new UnknownFunction(new Lagrange(2));
      /* dx is the the differentiation operator in direction 0. */
      Expr dx = new Derivative(0);
			
      /* create a discrete space on the mesh */
      TSFVectorSpace discreteSpace 
        = new SundanceVectorSpace(mesh, new Lagrange(2));

      /* create an expression for the initial guess. This will be reused as the
       * starting point for each newton step. Assume T(x)=x as an initial 
       * guess, and discretize it.
       */
      double tLeft = 1.0;
      double tRight = 2.0;
      Expr T0 = new DiscreteFunction(discreteSpace, tLeft + x*(tRight-tLeft));


      /* linearization of the weak PDE */
      Expr linearizedExpr = (dx*varT)*(dx*(pow(T0, 4) + 4.0*pow(T0, 3)*dT));
      Expr eqn = Integral(linearizedExpr);

      /* dirichlet BCs at both ends */
      EssentialBC bc = 
        EssentialBC(left, (T0+dT-tLeft)*varT) 
        && EssentialBC(right, (T0+dT-tRight)*varT) ;
			
      /* linear problem for the step du */
      StaticLinearProblem prob(mesh, eqn, bc, varT, dT); 

      /* create linear solver */
      TSFPreconditionerFactory prec = new ILUKPreconditionerFactory(1);
      TSFLinearSolver solver = new BICGSTABSolver(1.0e-14, 500);

      /* use Newton's method with backtracking to do the nonlinear solve */
			
      int maxIters = 100;
      int maxBacktracks = 10;
      double tol = 1.0e-12;
      double resid;
      double tResid = 1.0;
      int iter=0;

      TSFVector currentPoint;
      TSFVector oldPoint;
      T0.getVector(currentPoint);


      /* main newton loop */
      while (tResid > tol && iter < maxIters)
        {
          cerr << "newton step " << iter << endl;

          /* copy the initial point for this step. 
           * Each trial step will refer to this point. */
          oldPoint = currentPoint.deepCopy();

          /* compute the residual at the initial point */
          double initialResid = prob.getRHS().norm2();
          cerr << "initial residual = " << initialResid << endl;

          /* if the initial residual is small, we're done */
          if (initialResid < tol) 
            {
              tResid = 0.0;
              continue;
            }

          /* solve linearized problem for the newton step */
          Expr fullStep= prob.solve(solver);

          /* get the step in TSF form, so we can do vector operations on it */
          TSFVector stepVec;
          fullStep.getVector(stepVec);

          /* backtracking loop. First try the full Newton step. If the 
           * residual does not decrease, shorten the step size until
           * the residual decreases. */

          double step = 1.0;
          bool decrease = false;
          int backtrack = 0;
          while (!decrease)
            {
              backtrack++;
              if (backtrack > maxBacktracks) TSFError::raise("too many backtracks");

              /* move from the initial point by a fraction of the
                 Newton step. */
              currentPoint = oldPoint + step*stepVec;

              /* get the residual at the new point. This can be
               * obtained by evaluating the RHS of the linear problem
               * */
              T0.setVector(currentPoint);
              /* tell the problem that we need to recompute the RHS */
              prob.flushMatrixValues();
              /* compute the RHS and its norm */
              TSFVector residual = prob.getRHS();
              resid = residual.norm2();
              cerr << "residual norm = " << resid << endl;

              /* check for a decrease. */
              if (resid < initialResid) 
                {
                  decrease = true;
                }
              else
                {
                  decrease = false;
                  step = 0.5*step;
                }
            }

          iter++;
          /* if the step size is small enough, we're done */
          tResid = step*stepVec.norm2();
          cerr << "step resid " << tResid << endl;
        }

      /* check for convergence failure */
      if (iter >= maxIters)
        {
          TSFError::raise("backtracking newton failed to converge");
        }

      /* whew, we made it. */
      cerr << "backtracking newton converged after " << iter << " iterations "
           << endl;


      /* Now compare exact and computed solns. */
      Expr exactSoln = pow(pow(tLeft,4.0) 
                           + (pow(tRight,4.0)-pow(tLeft,4.0))*x, 0.25);
      Expr errorExpr = T0 - exactSoln;
			

      double errorNorm = errorExpr.norm(2, new GaussianQuadrature(6));
      double tolerance = 1.0e-4;
			
      Testing::passFailCheck(__FILE__, errorNorm, tolerance);
      Testing::timeStamp(__FILE__, __DATE__, __TIME__);

    }
  catch(exception& e)
    {
      TSFOut::println(e.what());
      Testing::crash(__FILE__);
      Testing::timeStamp(__FILE__, __DATE__, __TIME__);
    }
  Sundance::finalize();
}