eminisolve.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2009 Werner Scholz

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    magpar is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

static char const Id[] = "$Id: eminisolve.c 2795 2009-10-19 18:54:30Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "emini.h"
#include "field/field.h"
#include "util/util.h"
#include "llg/llg.h"
#include "io/magpario.h"

#ifdef ADDONS
#include "addons/addons.h"
#endif

#undef DEBUG
#ifdef DEBUG
static Vec dMdt=PETSC_NULL;
#endif

/* undef DAMP to use conventional gradient
   define DAMP to use damping term for gradient/line search
*/
#undef DAMP

/* TODO: problems with using Cartesian coordinates in micromagnetic model:
   |M|==1 is not preserved! M "runs away"! - need penalty term or such
   works great with nonlinear model
*/

static TAO_SOLVER tao;
static TAO_APPLICATION taoapp;
static PetscRandom rctx=PETSC_NULL;
static PetscReal magdist;
static PetscReal tol;
static int tao_max_its;

static Vec Msphere;
static Vec VHtotsphere;

#undef TESTTORQ
#ifdef TESTTORQ
static Vec dMdt=PETSC_NULL;
#endif

#ifdef DAMP
static Vec grad;
#endif

int TaoEvalEnergyGrad(TAO_APPLICATION taoapp,Vec X,double *f, Vec G,void *ptr)
{
  GridData *gdata = (GridData*)ptr;

  MagparFunctionInfoBegin;

  if (gdata->mode==53) {
    /* FIXME: try to avoid copying vectors here
       can we do gdata->M=X instead?
       same for gdata->VHtot=G ?
       search for all VecCopy and try to do something smarter
    */
    ierr = VecCopy(X,gdata->M);CHKERRQ(ierr);
  }
  else {
    ierr = Sphere2Cart(X,gdata->M);CHKERRQ(ierr);
  }

  ierr = Htot_Gradient(gdata);CHKERRQ(ierr);

#ifdef DAMP
  /* use damping term instead of field for gradient/line search */

  PetscReal *ta_M, *ta_Htot, *ta_grad;
  ierr = VecGetArray(gdata->M,&ta_M);CHKERRQ(ierr);
  ierr = VecGetArray(gdata->VHtot,&ta_Htot);CHKERRQ(ierr);
  ierr = VecGetArray(grad,&ta_grad);CHKERRQ(ierr);

  for (int i=0; i<gdata->ln_vert; i++) {
    PetscReal  *damp, mxh[ND];
    PetscReal  a1;

    /* skip non-magnetic vertices */
    if (gdata->propdat[NP*gdata->vertprop[i]+4]<=0.0) continue;

    damp=ta_grad+ND*i;

    douter(ND,ta_M+ND*i,ta_Htot+ND*i,mxh);
    douter(ND,ta_M+ND*i,mxh,damp);

    /* use small parameter to keep residual for TAO small */
    a1=0.001;
    my_dscal(ND,a1,damp,1);
  }

  ierr = VecRestoreArray(gdata->M,&ta_M);CHKERRQ(ierr);
  ierr = VecRestoreArray(gdata->VHtot,&ta_Htot);CHKERRQ(ierr);
  ierr = VecRestoreArray(grad,&ta_grad);CHKERRQ(ierr);
#endif

  ierr = Htot_Energy(gdata);CHKERRQ(ierr);
  *f=gdata->Etot;

/* Test gradient with these options:
  -tao_method tao_fd_test
  -tao_test_gradient
  -tao_test_display
*/

#ifdef DAMP
  if (gdata->mode==53) {
    ierr = VecCopy(grad,G);CHKERRQ(ierr);
  }
  else {
    ierr = Cart2SphereDiff(grad,X,G);CHKERRQ(ierr);
  }
#else
  if (gdata->mode==53) {
    ierr = VecCopy(gdata->VHtot,G);CHKERRQ(ierr);
  }
  else {
    ierr = Cart2SphereDiff(gdata->VHtot,X,G);CHKERRQ(ierr);
  }
#endif

  /* DEBUG
  PetscPrintf(PETSC_COMM_WORLD,"deb01 X\n");
  VecView(X,PETSC_VIEWER_STDOUT_WORLD);
  PetscPrintf(PETSC_COMM_WORLD,"deb01 VHtot\n");
  VecView(gdata->VHtot,PETSC_VIEWER_STDOUT_WORLD);
  PetscPrintf(PETSC_COMM_WORLD,"deb01 G\n");
  VecView(G,PETSC_VIEWER_STDOUT_WORLD);
  */

  MagparFunctionProfReturn(0);
}

int ConvTest(TAO_SOLVER tao,void *cctx)
{
  GridData *gdata = (GridData*)cctx;

  MagparFunctionInfoBegin;

  int iter;
  double ff,gnorm,cnorm,xdiff;
  TaoTerminateReason reason;
  ierr = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,&cnorm,&xdiff,&reason);CHKERRQ(ierr);
  if (iter>tao_max_its) {
    ierr = TaoSetTerminationReason(tao,TAO_DIVERGED_MAXITS);CHKERRQ(ierr);
    MagparFunctionProfReturn(0);
  }

#ifdef TESTTORQ
  /* test based on torque dM/dt - does not work well
     torque>0 in energy minimum!?
     however torque->0 is not (strict) energy minimum
  */
  if (gdata->mode!=53) {
    Vec sol;
    ierr = TaoAppGetSolutionVec(taoapp,&sol);CHKERRQ(ierr);
    assert(sol==Msphere);
    ierr = Sphere2Cart(sol,gdata->M);CHKERRQ(ierr);
  }

  ierr = calc_dMdt(PETSC_NULL,0.0,gdata->M,dMdt,gdata);CHKERRQ(ierr);
  gdata->vequil=RenormVec(dMdt,0.0,PETSC_MAX,ND);

  ierr = TaoGetTerminationReason(tao,&reason);CHKERRQ(ierr);

  PetscPrintf(PETSC_COMM_WORLD,"torque = %g  Etot = %g %g  reason = %i\n",
    gdata->vequil,gdata->Etot,gdata->Etot*gdata->escale,reason
  );

  static int equil=0;
  if (gdata->vequil < tol) {
    ierr = TaoSetTerminationReason(tao,TAO_CONVERGED_USER);CHKERRQ(ierr);
    /* count how long we have been in equilibrium in succession */
    equil++;
  }
  else {
    ierr = TaoSetTerminationReason(tao,TAO_CONTINUE_ITERATING);CHKERRQ(ierr);
    /* reset counter if we are not in equilibrium */
    equil=0;
  }

#else
  /* test based on change in energy */
  static PetscReal Etotold;
  static int equil;

  if (iter==0) {
    equil=0;
    Etotold=PETSC_MAX;
  }

  PetscLogDouble t_now;
  PetscGetTime(&t_now);
  static PetscLogDouble t_last=t_now;

  PetscPrintf(PETSC_COMM_WORLD,"%5i %4i %6.2f %22.16e %22.16e %22.16e %22.16e %22.16e",
    iter,
    equil,
    t_now-t_last,
    gdata->Etot,gdata->Etot*gdata->escale,
    Etotold,Etotold*gdata->escale,
    (Etotold-gdata->Etot)/PetscAbsReal(gdata->Etot)
  );
#ifdef DEBUG
  ierr = calc_dMdt(PETSC_NULL,0.0,gdata->M,dMdt,gdata);CHKERRQ(ierr);
  gdata->vequil=RenormVec(dMdt,0.0,PETSC_MAX,ND);

  PetscPrintf(PETSC_COMM_WORLD,"%22.16e",
    gdata->vequil
  );
#endif
  PetscPrintf(PETSC_COMM_WORLD,"\n");
  t_last=t_now;

  /* Etotold-gdata->Etot < 0 when energy increases
     This can happen around energy minimum, due to numerical errors.
     We allow it as long as |delta E|<tol.
  */
#define NEQUIL 5
  if (PetscAbs(Etotold-gdata->Etot)/PetscAbsReal(gdata->Etot)<tol) {
    if (equil>=NEQUIL) {
      ierr = TaoSetTerminationReason(tao,TAO_CONVERGED_USER);CHKERRQ(ierr);
      equil=0;
    }
    else {
      ierr = TaoSetTerminationReason(tao,TAO_CONTINUE_ITERATING);CHKERRQ(ierr);
      /* count how long we have been in equilibrium in succession */
      equil++;
    }
  }
  else {
    ierr = TaoSetTerminationReason(tao,TAO_CONTINUE_ITERATING);CHKERRQ(ierr);
    /* reset counter if we are not in equilibrium */
    equil=0;
  }
  Etotold=gdata->Etot;
#endif

#ifdef DEBUG
  ierr = WriteLog(gdata);CHKERRQ(ierr);
  gdata->inp=PetscAbsInt(gdata->inp);
  ierr = WriteSet(gdata);CHKERRQ(ierr);
#endif

  MagparFunctionProfReturn(0);
}


int myTSCreateEmini(GridData *gdata)
{
  MagparFunctionLogBegin;

  gdata->time = 0.0;

#ifdef DAMP
  PetscPrintf(PETSC_COMM_WORLD,"DAMP defined: Using damping term for gradient/line search\n");
  ierr = VecDuplicate(gdata->M,&grad);CHKERRQ(ierr);
#else
  PetscPrintf(PETSC_COMM_WORLD,"DAMP undefined: Using field for gradient/line search\n");
#endif

  /* set default for tao_lmm_vectors */
  PetscTruth flg;
  int tao_lmm_vectors;
  ierr = PetscOptionsGetInt(PETSC_NULL,"-tao_lmm_vectors",(PetscInt*)&tao_lmm_vectors,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    ierr = PetscOptionsSetValue("-tao_lmm_vectors","100");CHKERRQ(ierr);
  }

  char str[256];
  PetscOptionsGetString(PETSC_NULL,"-tao_method",str,256,&flg);
  if (flg!=PETSC_TRUE) {
    ierr = PetscSNPrintf(str,255,"tao_lmvm");CHKERRQ(ierr);
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -tao_method not found, using default value: %s\n",
      str
    );
  }

  ierr = TaoCreate(PETSC_COMM_WORLD,str,&tao);CHKERRQ(ierr);
  ierr = TaoApplicationCreate(PETSC_COMM_WORLD,&taoapp);CHKERRQ(ierr);
  if (gdata->mode==53) {
    PetscPrintf(PETSC_COMM_WORLD,"mode %i: Using 3D Cartesian coordinates\n",gdata->mode);

    PetscReal *ta_M;
    ierr = VecGetArray(gdata->M,&ta_M);CHKERRQ(ierr);
    for (int j=0;j<gdata->ln_vert;j++) {
      PetscReal ms;
      ms=gdata->propdat[NP*gdata->vertprop[j]+4];
      /* if the magnetic moment of the material is 0 */
      if (ms==0.0) {
        /* set magnetization to 0 */
        /* This is necessary, because we set Hnonlin=0 for Ms=0 or mu=1 / chi=0 in hnonlin.c */
        ta_M[ND*j+0] = 0.0;
        ta_M[ND*j+1] = 0.0;
        ta_M[ND*j+2] = 0.0;
      }
    }
    ierr = VecRestoreArray(gdata->M,&ta_M);CHKERRQ(ierr);

    ierr = TaoAppSetInitialSolutionVec(taoapp,gdata->M);CHKERRQ(ierr);
  }
  else {
    PetscPrintf(PETSC_COMM_WORLD,"Using spherical coordinates\n");
    ierr = VecCreate(PETSC_COMM_WORLD,&Msphere);CHKERRQ(ierr);
    ierr = VecSetSizes(Msphere,(ND-1)*gdata->ln_vert,(ND-1)*gdata->n_vert);CHKERRQ(ierr);
    ierr = VecSetFromOptions(Msphere);CHKERRQ(ierr);

    ierr = VecDuplicate(Msphere,&VHtotsphere);CHKERRQ(ierr);

    ierr = TaoAppSetInitialSolutionVec(taoapp,Msphere);CHKERRQ(ierr);
  }
  ierr = TaoAppSetObjectiveAndGradientRoutine(taoapp,TaoEvalEnergyGrad,(void *)gdata);CHKERRQ(ierr);
  ierr = TaoSetConvergenceTest(tao,ConvTest,(void *)gdata);CHKERRQ(ierr);

/*
  PetscPrintf(PETSC_COMM_WORLD,"Activating TAO Monitor\n");
  ierr = PetscOptionsSetValue("-tao_monitor","");CHKERRQ(ierr);
*/

  /* set defaults for line search options */
  ierr = PetscOptionsHasName(PETSC_NULL,"-tao_ls_ftol",&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    ierr = PetscOptionsSetValue("-tao_ls_ftol","1e-20");CHKERRQ(ierr);
  }
  ierr = PetscOptionsHasName(PETSC_NULL,"-tao_ls_rtol",&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    ierr = PetscOptionsSetValue("-tao_ls_rtol","1e-1");CHKERRQ(ierr);
  }
  ierr = PetscOptionsHasName(PETSC_NULL,"-tao_ls_gtol",&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    ierr = PetscOptionsSetValue("-tao_ls_gtol","0.99");CHKERRQ(ierr);
  }

  ierr = PetscOptionsGetInt(PETSC_NULL,"-tao_max_its",(PetscInt*)&tao_max_its,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    tao_max_its=500;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -tao_max_its not found, using default value: %i\n",
      tao_max_its
    );
  }

  /* Check for TAO command line options (listed here for consistency check with allopt.txt) */
  /* evaluates
      PetscOptions "-tao_fatol"
      PetscOptions "-tao_frtol"
  */
/* defined in tao-1.9/src/petsctao/include/taopetsc.h - may not be supported in future releases
  ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);
*/
  ierr = TaoSetOptions(taoapp,tao);CHKERRQ(ierr);

  /*
    disturb magnetization distribution with small random vector;
    larger distortions require many more interations in TaoSolve
    and do not help, because the solution is less accurate
    (relative tolerance!)
  */
  ierr = PetscOptionsGetReal(PETSC_NULL,"-magdist", &magdist,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    magdist=0.0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -magdist not found, using default value: %g\n",
      magdist
    );
  }

  if (magdist<0.0) {
    PetscPrintf(PETSC_COMM_WORLD,"magdist = %g < 0.0: Distorting M by tilting in direction MxMxH\n",magdist);
  }
  else if (magdist>0.0) {
    PetscPrintf(PETSC_COMM_WORLD,"magdist = %g > 0.0: Distorting M randomly\n",magdist);
  }
  else {
    PetscPrintf(PETSC_COMM_WORLD,"magdist==0.0: Not distorting M\n");
  }

  ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
  ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);

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

  ierr = PetscOptionsGetReal(PETSC_NULL,"-tol",&tol,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    tol=1e-7;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -tol not found, using default value: %g\n",
      tol
    );
  }
  PetscPrintf(PETSC_COMM_WORLD,"tol (equilibrium): %g\n",tol);

#ifdef TESTTORQ
  ierr = VecDuplicate(gdata->M,&dMdt);
#endif
#ifdef DEBUG
  ierr = VecDuplicate(gdata->M,&dMdt);
#endif

  MagparFunctionLogReturn(0);
}


int EminiSolve(GridData *gdata)
{
  MagparFunctionInfoBegin;

  if (magdist<0.0) {
    /* distort by tilting magnetization in direction of MxMxH */
    PetscReal *ta_M, *ta_Htot;
    VecGetArray(gdata->M,&ta_M);
    VecGetArray(gdata->VHtot,&ta_Htot);

    for (int i=0; i<gdata->ln_vert; i++) {
      PetscReal  damp[ND], mxh[ND];

      douter(ND,ta_M+ND*i,ta_Htot+ND*i,mxh);
      douter(ND,ta_M+ND*i,mxh,damp);

      /* TODO: can lead to Inf/NaNs if "damping constant" is too large */
      my_daxpy(ND,magdist,damp,1,ta_M+ND*i,1);
    }

    VecRestoreArray(gdata->M,&ta_M);
    VecRestoreArray(gdata->VHtot,&ta_Htot);
  }
  else if (magdist>0.0) {
    /* randomly distort M */
    ierr = DistortVec(rctx,gdata->M,magdist);CHKERRQ(ierr);

    /* renormalize M */
    if (gdata->mode!=53) {
      PetscReal devNorm;
      devNorm=RenormVec(gdata->M,1.0,D_EPS,ND);
    }
  }

  PetscPrintf(MPI_COMM_WORLD,"Starting energy minimization with TAO...\n");
  PetscPrintf(PETSC_COMM_WORLD,"iter  equil  CPUtstep (s)   Etot (1)            (J/m^3)                Etotold  (1)             (J/m^3)               tol               torq\n");
  PetscPrintf(PETSC_COMM_WORLD,"--------------------------------------------------------------------------------------------------------------------------------------------------------------\n");

  if (gdata->mode!=53) {
    ierr = Cart2Sphere(gdata->M,Msphere);CHKERRQ(ierr);
  }

  ierr = TaoSolveApplication(taoapp,tao);CHKERRQ(ierr);

#ifdef ADDONS

#ifdef DEBUG
  /* recalculate fields and energies */
  ierr = Htot(gdata);CHKERRQ(ierr);
  ierr = Htot_Energy(gdata);CHKERRQ(ierr);

  gdata->inp=PetscAbsInt(gdata->inp);
  ierr = WriteLog(gdata);CHKERRQ(ierr);
  ierr = WriteSet(gdata);CHKERRQ(ierr);
#endif

  PetscTruth flg;
  int NSMOOTHSTEPS;
  ierr = PetscOptionsGetInt(PETSC_NULL,"-emini_mag_smooth_steps",(PetscInt*)&NSMOOTHSTEPS,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    NSMOOTHSTEPS=1;
  }

  if (NSMOOTHSTEPS>0) {
    PetscPrintf(MPI_COMM_WORLD,"Smoothening magnetization distribution\n");

    /* smoothen magnetization distribution */
    Vec Vtmp;
    ierr = VecDuplicate(gdata->M,&Vtmp);CHKERRQ(ierr);

    for (int i=0;i<NSMOOTHSTEPS;i++) {
      /* perform smoothening (double) step */
      ierr = FieldSmoothApply(gdata,gdata->M,Vtmp);
      ierr = FieldSmoothApply(gdata,Vtmp,gdata->M);

      /* renormalize M (distorted due to smoothening) */
      PetscReal devNorm;
      devNorm=RenormVec(gdata->M,1.0,0.0,ND);
    }

    ierr = VecDestroy(Vtmp);CHKERRQ(ierr);
  }

  if (gdata->mode!=53) {
    ierr = Cart2Sphere(gdata->M,Msphere);CHKERRQ(ierr);
  }

  if (NSMOOTHSTEPS>0) {
    PetscPrintf(MPI_COMM_WORLD,"Restarting energy minimization with TAO after smoothening...\n");
    ierr = TaoSolveApplication(taoapp,tao);CHKERRQ(ierr);
  }
#endif

  if (gdata->mode!=53) {
    Vec sol;
    ierr = TaoAppGetSolutionVec(taoapp,&sol);CHKERRQ(ierr);
    assert(sol==Msphere);
  }
  /* store this info in log file!? */
  int iter;
  double ff,gnorm,cnorm,xdiff;
  TaoTerminateReason reason;
  ierr = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,&cnorm,&xdiff,&reason);CHKERRQ(ierr);
  if (reason <= 0 ){
    PetscPrintf(MPI_COMM_WORLD,"TAO solver did NOT converge - ");
    /* continue even if TAO fails
       SETERRQ(PETSC_ERR_LIB,"TAO failed! Exiting!\n");
    */
    gdata->equil=0;
  }
  else {
    PetscPrintf(MPI_COMM_WORLD,"TAO solver converged - ");
    gdata->equil++;
  }
  PetscPrintf(MPI_COMM_WORLD,
    "TerminationReason: %i\n"
    "it: %i, f: %4.2e, res: %4.2e, cnorm: %4.2e, step %4.2e\n",
    reason,
    iter,ff,gnorm,cnorm,xdiff
  );

  ierr = TaoView(tao);CHKERRQ(ierr);

  /* otherwise Htot=0 */
  ierr = Htot(gdata);CHKERRQ(ierr);
  ierr = Htot_Energy(gdata);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}