/*$Id: ginzland3.c,v 1.1 2002/05/13 19:32:57 benson Exp $*/

/* Program usage: mpirun -np <proc> ginzland3 [-help] [all TAO options] */

/*
  Include "tao.h" so we can use TAO solvers.
  petscda.h for distributed array
  ad_deriv.h for AD gradient
*/

#include "petscda.h"
#include "tao.h"
#include "taodaapplication.h"


static char  help[] = "";

/*
   User-defined application context - contains data needed by the
   application-provided call-back routines.
*/


/*  
    This structure is used only when an ADIC generated gradient is used.
    An InactiveDouble type is a double 
*/
typedef struct {
  PetscScalar x,y,vpotx,vpoty;
} Field;

typedef struct {
  int                 mx,my;
  InactiveDouble      hx, hy;        /* increment size in both directions */
  InactiveDouble      twopivornum; 
  InactiveDouble      tkappa;
} ADFGCtx;

#ifdef ad_GRAD_MAX
#undef ad_GRAD_MAX
#endif
#define ad_GRAD_MAX 16

int ad_GinzLandLocalFunction(int[2], DERIV_TYPE[16], DERIV_TYPE*, void*);
/*
int DAAppSetADElementFunctionGradient2(TAO_APPLICATION, 
				       int (*)(int[2],DERIV_TYPE[16],DERIV_TYPE*,void*), 
				       int, void *);
*/
typedef struct {
  int         mx,my;
  int         vornum;
  double      tkappa;
  double      hx, hy;        /* increment size in both directions */
  ADFGCtx     fgctx;         /* Used only when an ADIC generated gradient is used */
} AppCtx;

/* User-defined routines */
static int AppCtxInitialize(void *);
static int AppInitialSolution(TAO_APPLICATION myapp, DA da, Vec X, AppCtx *ctx);
static int MyGridMonitorBefore(TAO_APPLICATION, DA, int, void *);

#undef __FUNCT__
#define __FUNCT__ "main"
int main( int argc, char **argv ) {

  int                  info;         /* used to check for functions returning nonzeros */
  int                  iter,nlevels;
  int                  Nx,Ny;
  double               ff,gnorm;
  DA                   DAarray[20];
  PetscTruth           flg, PreLoad = PETSC_TRUE;                             /* flags */
  TaoMethod            method = "tao_nls";                      /* minimization method */
  TaoTerminateReason   reason;
  TAO_SOLVER           tao;                               /* TAO_SOLVER solver context */
  TAO_APPLICATION      GinzburghLandau;                       /* The PETSc application */
  AppCtx               user;                              /* user-defined work context */
  Vec x;
  double zero=0.0;

  /* Initialize TAO */
  PetscInitialize(&argc, &argv, (char *)0, help);
  TaoInitialize(&argc, &argv, (char *)0, help);

  PreLoadBegin(PreLoad,"Solve");

  info = AppCtxInitialize((void*)&user); CHKERRQ(info);

  nlevels=1;
  info = PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,&flg); CHKERRQ(info);

  if (PreLoadIt == 0) {
    nlevels = 1; user.mx = 12; user.my = 12;
  }

  PetscPrintf(MPI_COMM_WORLD,"\n---- Ginzburg - Landau Superconductivity -----\n\n");

  /* Let PETSc determine the vector distribution */
  Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

  /* Create distributed array (DA) to manage parallel grid and vectors  */
  info = DACreate2d(PETSC_COMM_WORLD,DA_XYPERIODIC, DA_STENCIL_BOX, user.mx,
                    user.my,Nx,Ny,4,1,PETSC_NULL,PETSC_NULL,&DAarray[0]); CHKERRQ(info);
  for (iter=1;iter<nlevels;iter++){
    info = DARefine(DAarray[iter-1],PETSC_COMM_WORLD,&DAarray[iter]); CHKERRQ(info);
  }

  /* Create TAO solver and set desired solution method */
  info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);
  info = TaoApplicationCreate(PETSC_COMM_WORLD, &GinzburghLandau); CHKERRQ(info);
  info = TaoAppSetDAApp(GinzburghLandau,DAarray,nlevels); CHKERRQ(info);
  info = DAAppSetBeforeMonitor(GinzburghLandau,MyGridMonitorBefore,(void*)&user); CHKERRQ(info);

  info = TaoAppSetRelativeTolerance(GinzburghLandau,1.0e-8); CHKERRQ(info);
  info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
  info = TaoSetGradientTolerances(tao,0,0,0); CHKERRQ(info);
  info = DAAppSetADElementFunctionGradient2(GinzburghLandau,ad_GinzLandLocalFunction,150,(void *)&user.fgctx); CHKERRQ(info);
  info = DAAppSetHessianMat(GinzburghLandau); CHKERRQ(info);
  info = TaoAppSetHessianRoutine(GinzburghLandau,TaoAppDefaultComputeHessianColor,(void*)GinzburghLandau);CHKERRQ(info);

  /* Check for any tao command line options */
  info = DAAppGetSolution(GinzburghLandau,0,&x); CHKERRQ(info);
  info = VecSet(x, zero); CHKERRQ(info);
  info = AppInitialSolution(GinzburghLandau,DAarray[0], x, &user); CHKERRQ(info);
  //  info=VecView(x,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
  //  info = DAAppSetInitialSolution(GinzburghLandau, x); CHKERRQ(info);

  info = TaoSetOptions(GinzburghLandau,tao); CHKERRQ(info);

  /* SOLVE THE APPLICATION */
  info = TaoDAAppSolve(GinzburghLandau,tao);  CHKERRQ(info);

  /* Get information on termination */
  info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
  if (reason <= 0 ){
    PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
    PetscPrintf(MPI_COMM_WORLD," Iterations: %d,  Function Value: %4.2e, Residual: %4.2e \n",iter,ff,gnorm);
  }

  info = PetscOptionsHasName(PETSC_NULL,"-view_sol",&flg); CHKERRQ(info);
  if (flg){
    info = DAAppGetSolution(GinzburghLandau,nlevels-1,&x); CHKERRQ(info);
    info = VecView(x,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
  }

  /*  To View TAO solver information */
  info = TaoView(tao); CHKERRQ(info);

  /* Free TAO data structures */
  info = TaoDestroy(tao); CHKERRQ(info);
  info = TaoAppDestroy(GinzburghLandau); CHKERRQ(info);

  /* Free PETSc data structures */
  for (iter=0;iter<nlevels;iter++){
    info = DADestroy(DAarray[iter]); CHKERRQ(info);
  }

  PreLoadEnd();

  /* Finalize TAO */
  TaoFinalize();
  PetscFinalize();

  return 0;
} /* main */

static int AppInitialSolution(TAO_APPLICATION myapp, DA da, Vec X, AppCtx *ctx) {

  int xs, xm, mx, ys, ym, my;
  double tkappa, hx, hy, sqn, ppi,bave, sqrtv, vornum;
  double xpt,ypt,t1;
  Field **x;
  int i,j, info;

  vornum=ctx->vornum;
  tkappa=ctx->tkappa;

  info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
                   PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);
  info = DAGetCorners(da, &xs, &ys, PETSC_NULL, &xm, &ym, PETSC_NULL); CHKERRQ(info);

  info = MyGridMonitorBefore(myapp,da,0,(void*)ctx); CHKERRQ(info);

  tkappa=ctx->tkappa; vornum=ctx->vornum;
  hx=ctx->hx; hy=ctx->hy;
  sqn=mx*my; ppi=4.0*atan(1.0);
  bave=2.0*(ppi)*(vornum)*(tkappa)/(sqn*hx*hy);
  sqrtv=sqrt(vornum*1.0)*(ppi);

  info = DAVecGetArray(da, X, (void**)&x); CHKERRQ(info);
  for (j = ys; j < ys+ym; j++) {
    ypt=j*hy;
    for (i = xs; i < xs+xm; i++) {
      xpt=i*hx;
      t1= sin(sqrtv*xpt/6.0) * sin(sqrtv*ypt/(6.0*sqrt(3.0))) ;
      x[j][i].x = 1.0 - t1*t1;
      x[j][i].y = 0.0;
    }
  }

  for (j = ys; j < ys+ym; j++) {
    for (i = xs; i < xs+xm; i++) {
      xpt=i*hx;
      x[j][i].vpotx = 0.0;
      x[j][i].vpoty = bave*xpt/tkappa;
    }
  }
  info = DAVecRestoreArray(da, X, (void**)&x); CHKERRQ(info);

  return 0;
}

/*----- The following two routines
  MyGridMonitorBefore    MyGridMonitorAfter
  help diplay info of iterations at every grid level 
*/

#undef __FUNCT__
#define __FUNCT__ "MyGridMonitorBefore"
static int MyGridMonitorBefore(TAO_APPLICATION myapp, DA da, int level, void *ctx) {

  AppCtx *user = (AppCtx*)ctx;
  int info,mx,my;
  double vornum;

  info = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
		   PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(info);

  vornum=user->vornum;
  user->hx=sqrt(vornum/2.0)*3.0/(mx);
  user->hy=sqrt(vornum/2.0)*3.0*sqrt(3.0)/(my);

  user->mx = mx;
  user->my = my;

  user->fgctx.hx          = user->hx;
  user->fgctx.hy          = user->hy;
  user->fgctx.mx          = user->mx;
  user->fgctx.my          = user->my;
  user->fgctx.tkappa      = user->tkappa;
  user->fgctx.twopivornum = 2.0*(4.0*atan(1.0))*vornum;
  { 
    KSP ksp=0;
    PC   pc=0;
    info = TaoAppGetKSP(myapp,&ksp);CHKERRQ(info);
    if (ksp){
      int maxits;
      info = KSPGetPC(ksp,&pc);CHKERRQ(info);
      if (mx <= 40){
	info = KSPGetTolerances(ksp,0,0,0,&maxits); CHKERRQ(info);
	info = KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxits*2);CHKERRQ(info);
      } else {
	info = KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,mx*10);CHKERRQ(info);
      }
      if (level <0){
	info = PCSetType(pc,PCJACOBI);CHKERRQ(info);
	info = KSPSetType(ksp,KSPGMRES); CHKERRQ(info);
      } else {
	info = PCSetType(pc,PCJACOBI);CHKERRQ(info);
	info = KSPSetType(ksp,KSPGMRES); CHKERRQ(info);
      }
    }
  }

  return 0;
}


#undef __FUNCT__
#define __FUNCT__ "AppCtxInitialize"
/*
  AppCtxInitialize - Sets initial values for the application context parameters

  Input:
    ptr - void user-defined application context

  Output:
    ptr - user-defined application context with the default or user-provided
             parameters
*/
static int AppCtxInitialize(void *ptr) {

  AppCtx *user = (AppCtx*)ptr;
  PetscTruth flg;
  int info;

  /* Specify default parameters */
  user->mx = user->my = 12;
  user->vornum=8;
  user->tkappa=5.0;

  /* Check for command line arguments that override defaults */
  info = PetscOptionsGetInt(TAO_NULL, "-mx", &user->mx, &flg); CHKERRQ(info);
  info = PetscOptionsGetInt(TAO_NULL, "-my", &user->my, &flg); CHKERRQ(info);
  info = PetscOptionsGetInt(TAO_NULL, "-par", &user->vornum, &flg); CHKERRQ(info);
  info = PetscOptionsGetReal(TAO_NULL, "-kappa", &user->tkappa, &flg); CHKERRQ(info);

  user->hx=sqrt(user->vornum/2.0)*3.0/(user->mx);
  user->hy=sqrt(user->vornum/2.0)*3.0*sqrt(3.0)/(user->my);

  user->fgctx.hx          = user->hx;
  user->fgctx.hy          = user->hy;
  user->fgctx.mx          = user->mx;
  user->fgctx.my          = user->my;
  user->fgctx.tkappa      = user->tkappa;
  user->fgctx.twopivornum = 2.0*(4.0*atan(1.0))*(user->vornum);

  PetscLogFlops(8);

  return 0;

} /* AppCtxInitialize */


