#include "Sundance.h"


/** \example inlinePoissonBoltzmann1D.cpp
 * Solve the Poisson-Boltzmann equation \f$\nabla^2 u = e^-u$ on the unit
 * line with boundary conditions:
 * Left: Natural, du/dx=0
 * Right: Dirichlet u = 2 log(cosh(1/sqrt(2)))
 *
 * The solution is 2 log(cosh(x/sqrt(2))).
 *
 * The problem is nonlinear, so we use Newton's method to iterate 
 * towards a solution. In this example, Newton's method is 
 * coded inline; we don't use any sort
 * of nonlinear solver object. 
 *
 */

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

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

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

      /* create a cell set representing the right boundary */
      CellSet boundary = new BoundaryCellSet();
      CellSet right = boundary.subset( x == L );

      /* create a discrete space on the mesh */
      TSFVectorSpace discreteSpace 
        = new SundanceVectorSpace(mesh, new Lagrange(1));

      /* create an expression for the initial guess. This will be reused as the
       * starting point for each newton step. Assume u(x)=x as an initial 
       * guess, and discretize it.
       */
      Expr u0 = new DiscreteFunction(discreteSpace, x, "u0");

      /* create symbolic objects for test and unknown functions. At each newton
       * step we will solve a linearized equation for a step du, so our
       * unknown is du. */
      Expr v = new TestFunction(new Lagrange(1), "v");
      Expr du = new UnknownFunction(new Lagrange(1), "du");

      /* create a differential operator representing the x-derivative. */
      Expr dx = new Derivative(0);
      
      /* linearized weak equation for the step du */
      Expr newton = (dx*(u0+du))*(dx*v) + v*exp(-u0) - v*exp(-u0)*du;
      Expr eqn = Integral(newton);

      /* Dirichlet boundary condition */
      double uBC = 2.0*log(cosh(L/sqrt(2.0)));          
      EssentialBC bc = EssentialBC(right, (u0+du-uBC)*v) ;
      
      /* linear problem for the step du */
      StaticLinearProblem prob(mesh, eqn, bc, v, du); 

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

      
      int maxIters = 10;
      double tol = 1.0e-12;
      double resid = 1;
      int iter=0;
      while (resid > tol && iter < maxIters)
        {
          TSFOut::printf("Newton step# %d", iter);

          /* solve for the step and compute its norm */
          Expr duSoln = prob.solve(solver);
          resid = duSoln.norm(2);

          TSFOut::printf("resid = %g\n", resid);
          iter++;

          /* add the full step to the previous point. Store the result
          * in the expression u0, so that its value in the linear
          * problem is updated. */
          u0 = new DiscreteFunction(discreteSpace, u0+duSoln);

          /* write the current iterate to disk */
          char fName[100];
          sprintf(fName, "u%d.dat.%d", iter, MPIComm::world().getRank());
          
          FieldWriter writer = new MatlabWriter(fName);
          writer.writeField(u0);

          /* discard the matrix values, preserving the matrix structure. */
          prob.flushMatrixValues();
        }

      
      
      // compare to exact solution
      Expr exactSoln = 2.0*log(cosh(x/sqrt(2.0)));

      // compute the norm of the error
      double errorNorm = (exactSoln - u0).norm(2);
      double tolerance = 1.0e-4;
      TSFOut::printf("error = %g\n", errorNorm);

      Testing::passFailCheck(__FILE__, errorNorm, tolerance);
    }
  catch(exception& e)
    {
      Sundance::handleError(e, __FILE__);
    }
  Sundance::finalize();
}