checkiterationllg.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: checkiterationllg.c 2681 2009-07-31 04:30:53Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

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

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

static int         doinit=1;

static PetscReal   tol;
static PetscReal   renormtol=-1.0;
static int         condinp_equil;
static PetscReal   condinp_j;
static PetscReal   condinp_equil_j;
static PetscReal   condinp_t;
static PetscReal   *condinp_user=PETSC_NULL;
static int         userinptot = 0;

/* TODO: resolve this extern dependency */
extern KSP ksp_Au1;
extern KSP ksp_Au2;
extern realtype abstol, reltol; 
int adjtoldem(GridData *gdata, double fac)
{
  MagparFunctionInfoBegin;

  if (gdata->VHdem==PETSC_NULL) {
    MagparFunctionInfoReturn(0);
  }

  PetscReal rtol, atol, dtol;
  int maxits;

  ierr = KSPGetTolerances(ksp_Au1,&rtol,&atol,&dtol,(PetscInt*)&maxits);CHKERRQ(ierr);

  if (rtol < 1e-10) {
    MagparFunctionInfoReturn(0);
  }

  atol *= fac;
  rtol *= fac;
  ierr = KSPSetTolerances(ksp_Au1,rtol,atol,dtol,maxits);CHKERRQ(ierr);

  PetscPrintf(PETSC_COMM_WORLD,
    "at t=%g ns: ",gdata->time*gdata->tscale*1e9
  );

  PetscPrintf(PETSC_COMM_WORLD,
    "%s KSP reltol to %g , abstol to %g\n",
    fac>1.0 ? "Increased" : "Reduced",
    rtol,atol
  );

  ierr = KSPGetTolerances(ksp_Au2,&rtol,&atol,&dtol,(PetscInt*)&maxits);CHKERRQ(ierr);
  atol *= fac;
  rtol *= fac;
  ierr = KSPSetTolerances(ksp_Au2,rtol,atol,dtol,maxits);CHKERRQ(ierr);

  MagparFunctionInfoReturn(0);
}


int adjtol(GridData *gdata, double fac)
{
  MagparFunctionInfoBegin;

  if (reltol < 1e-10) {
    MagparFunctionInfoReturn(0);
  }

  PetscPrintf(PETSC_COMM_WORLD,
    "at t=%g ns: ",gdata->time*gdata->tscale*1e9
  );

  adjtoldem(gdata,fac);

  reltol *= fac;
  abstol *= fac;

  extern void *cvode_mem;
  long int nstep;
  int flag;
  flag=CVodeGetNumSteps(cvode_mem,&nstep);
  PetscPrintf(PETSC_COMM_WORLD,
    "%s PVODE reltol to %g , abstol to %g ; nstep since last reinit: %i\n",
    fac>1.0 ? "Increased" : "Reduced",
    reltol,abstol,nstep
  );

  ierr = PVodeReInit(gdata);CHKERRQ(ierr);

  MagparFunctionInfoReturn(0);
}


int EquilCheck(Vec M,PetscReal time,PetscReal *torque)
{
  static PetscReal tlast=time-0.00001;
  static Vec Mold=PETSC_NULL;

  MagparFunctionInfoBegin;

  /* initialize upon first call */
  if (Mold==PETSC_NULL) {
    ierr = VecDuplicate(M,&Mold);CHKERRQ(ierr);
    ierr = VecZeroEntries(Mold);CHKERRQ(ierr);

    /* copy M */
    ierr = VecCopy(M,Mold);CHKERRQ(ierr);
    /* return fake value since we have not done any timesteps yet
       just larger than the tolerance to keep the simulation running
    */
    *torque=tol*1.01;

    MagparFunctionInfoReturn(0);
  }

  ierr = VecAXPY(Mold,-1.0,M);CHKERRQ(ierr);

  /* calculate norm of delta M */
  *torque=RenormVec(Mold,0.0,PETSC_MAX,ND)/(time-tlast);

  ierr = VecCopy(M,Mold);CHKERRQ(ierr);

  tlast=time;

  MagparFunctionProfReturn(0);
}


int CheckIterationLLG_Init(GridData *gdata)
{
  PetscTruth flg;

  MagparFunctionLogBegin;
  char fn[256];
  char line[256];
  int rank;
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

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

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

  if (abstol > tol/10.0) {
    PetscPrintf(PETSC_COMM_WORLD,
      "Warning: abstol=%g should be < (tol=%g)/10.0!\n",
      abstol,
      tol
    );
  }

  ierr = PetscOptionsGetInt(PETSC_NULL,"-condinp_equil", (PetscInt*)&condinp_equil, &flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    condinp_equil=1;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -condinp_equil not found, using default value: %i\n",
      condinp_equil
    );
  }
  ierr = PetscOptionsGetReal(PETSC_NULL,"-condinp_j",&condinp_j,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    condinp_j=0.1;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -condinp_j not found, using default value: %g\n",
      condinp_j
    );
  }
  ierr = PetscOptionsGetReal(PETSC_NULL,"-condinp_equil_j",&condinp_equil_j,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    condinp_equil_j=0.0;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -condinp_equil_j not found, using default value: %g\n",
      condinp_equil_j
    );
  }
  ierr = PetscOptionsGetReal(PETSC_NULL,"-condinp_t",&condinp_t,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    condinp_t=PETSC_MAX;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -condinp_t not found, using default value: %g\n",
      condinp_t
    );
  }

  /* user times in file name specified by -condinp_user */
  ierr = PetscOptionsGetString(PETSC_NULL,"-condinp_file_t_ns",fn,sizeof(fn),&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -condinp_file_t_ns not found, user inp output times NOT activated\n"
    );
    userinptot = 0;
  }
  else {
    if (!rank) {
      FILE *fd=NULL;
      PetscPrintf(PETSC_COMM_WORLD,"Opening condinp_file_t_ns file %s\n",fn);
      ierr = PetscFOpen(PETSC_COMM_WORLD,fn,"r",&fd);CHKERRQ(ierr);
      if (!fd){
        SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",fn);
      }
      read_one_line(fd,line,sizeof(line),NULL);
      if (sscanf(line,"%i",&userinptot) != 1) {
        SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"could not read number of user inp values\n");
      }
      /* alloc memory to store user time values */
      ierr = PetscMalloc(userinptot*sizeof(PetscReal),&condinp_user);CHKERRQ(ierr);

      for (int i=0;i<userinptot;i++) {
        read_one_line(fd,line,sizeof(line),NULL);
        if (sscanf(line,"%lf",condinp_user + i) != 1) {
        SETERRQ1(PETSC_ERR_FILE_UNEXPECTED,"could not read user inp time value %i\n",i);
        }
        condinp_user[i] /= 1e9*gdata->tscale;
      }

      PetscPrintf(PETSC_COMM_WORLD,"Number of user inp output times read in: %i\n",userinptot);
      ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);
    }
    /* broadcast condinp_user data to all processors */
    ierr = MPI_Bcast(&userinptot,1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
    if (rank) {
      ierr = PetscMalloc(userinptot*sizeof(PetscReal),&condinp_user);CHKERRQ(ierr);
    }
    ierr = MPI_Bcast(condinp_user,userinptot,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  }
  condinp_t /= 1e9*gdata->tscale;

  MagparFunctionLogReturn(0);
}


int CheckIterationLLG(GridData *gdata)
{
  static PetscReal Etotold=PETSC_MAX;
  static PetscReal hextold=PETSC_MAX;
  static PetscReal condinp_last_t=PETSC_MIN;
  static PetscReal condinp_last_j=PETSC_MIN;
  static PetscReal condinp_last_equil_j=PETSC_MIN;

  const PetscReal  tolscale=3.0;

  static int nsteps=-1;
  static int csteps;

  static PetscReal Edemold=gdata->Edem;
  static PetscReal dEdem1=0;
  static PetscReal dEdem2=0;

  static PetscReal timelast=gdata->time;

  static int userinpindex = 0;

#undef CHKDMDT
#ifdef CHKDMDT
  static Vec dMdt;
#endif

  MagparFunctionInfoBegin;

  /* initialize if necessary */
  if (doinit) {
    ierr = CheckIterationLLG_Init(gdata);CHKERRQ(ierr);
#ifdef CHKDMDT
    VecDuplicate(gdata->M,&dMdt);
#endif
    doinit=0;
  }

#ifdef CHKDMDT
  /* This method could be more expensive than the old method because
     we recalculate all the fields when we call calc_dMdt!
     However, we would get access the old dMdt - then it would be cheaper.
  */
  PetscReal veq;
  calc_dMdt(PETSC_NULL,0.0,gdata->M,dMdt,gdata);
  veq=RenormVec(dMdt,0.0,PETSC_MAX,ND);
  gdata->equil=veq;

#if 1
  ierr = EquilCheck(gdata->M,gdata->time,&gdata->vequil);CHKERRQ(ierr);
#endif

  PetscPrintf(PETSC_COMM_WORLD,"torque = %g  %g  %g\n",veq,gdata->vequil,veq/gdata->vequil);
#else
  ierr = EquilCheck(gdata->M,gdata->time,&gdata->vequil);CHKERRQ(ierr);
#endif

  /* We could also use differential dMdt calculated in calc_dMdt through RHSfunction (see there).
     However, there is more noise than if we take the finite difference over one full time step.
  */

  if (gdata->vequil < tol) {
    /* count how long we have been in equilibrium in succession */
    gdata->equil++;
  }
  else {
    /* reset counter if we are not in equilibrium */
    gdata->equil=0;
  }

  /* check if total energy increased (which it should not)
     print warning only if external field has not changed
  */
  PetscReal hext;
  hext=Hexternal_hext();
  if (gdata->Etot > Etotold &&
      hext == hextold &&
      ((gdata->Etot-Etotold)/Etotold)>1e-8) {

    PetscPrintf(MPI_COMM_WORLD,
      "Warning (t=%g ns): Etot increased by %g (%g*Etot) from %g to %g J/m^3\n",
      gdata->time*gdata->tscale*1e9,
      (gdata->Etot-Etotold)*gdata->escale,
      (gdata->Etot-Etotold)/Etotold,
      Etotold*gdata->escale,
      gdata->Etot*gdata->escale
    );

/*
    ierr = adjtol(gdata,1.0/tolscale);CHKERRQ(ierr);
*/
  }
  Etotold=gdata->Etot;
  hextold=hext;

  PetscReal valMpHext;
  ierr = MpHext(gdata->M,&valMpHext);CHKERRQ(ierr);

  /* check if a set of output data should be written */
  if (PetscAbsReal(valMpHext-condinp_last_j) > condinp_j) {
    gdata->inp=PetscAbsInt(gdata->inp);

    condinp_last_t=gdata->time;
    condinp_last_j=valMpHext;
  }
  else if (PetscAbsReal(gdata->time-condinp_last_t) > condinp_t) {
    gdata->inp=PetscAbsInt(gdata->inp);

    condinp_last_t=gdata->time;
    condinp_last_j=valMpHext;
  }
  else if ( userinpindex < userinptot && gdata->time >= condinp_user[userinpindex]) {
      gdata->inp=PetscAbsInt(gdata->inp);

      condinp_last_t=gdata->time;
      condinp_last_j=valMpHext;
      userinpindex++;
  }
  else if (condinp_equil==1 && gdata->equil > 0) {
    /* check for condinp_equil_j (defaults to 0.0) */
    if (PetscAbsReal(valMpHext-condinp_last_equil_j) > condinp_equil_j) {
      gdata->inp=PetscAbsInt(gdata->inp);

      condinp_last_t=gdata->time;
      condinp_last_j=valMpHext;
      condinp_last_equil_j=valMpHext;
    }
  }

  if (nsteps<0) {
    PetscTruth flg;
    ierr = PetscOptionsGetInt(PETSC_NULL,"-ts_logsteps",(PetscInt*)&nsteps,&flg);CHKERRQ(ierr);
    if (flg!=PETSC_TRUE) {
      nsteps=1;
      PetscPrintf(PETSC_COMM_WORLD,
        "Option -ts_logsteps not found, using default value: %i\n",
        nsteps
      );
    }
    if (nsteps<1) {
      SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Option set: -ts_logsteps %i, but has to be >0\n",nsteps);
    }
    csteps=0;
  }

  if (csteps % nsteps == 0 || gdata->equil>0) {
    ierr = WriteLog(gdata);CHKERRQ(ierr);
    ierr = WriteSet(gdata);CHKERRQ(ierr);
    csteps=0;
  }
  csteps++;

#ifdef ADDONS
  /* move magnetization if requested */
  ierr = MovingField(gdata);CHKERRQ(ierr);
#endif

  extern void *cvode_mem;
  long int nstep;
  int flag;
  flag=CVodeGetNumSteps(cvode_mem,&nstep);

  /* update external field if system is in equilibrium */
  if (gdata->equil>0) {
    ierr = Hext_ho_hstep(gdata);CHKERRQ(ierr);

    /* check if we are still in equilibrium */
    if (gdata->equil>0) {
      /* that's because the field has not been changed */
      PetscPrintf(PETSC_COMM_WORLD,"gdata->equil %i >= 1\n",gdata->equil);
      gdata->mode=99;
    }
    /* restart PVODE because the external field is changed abruptly */
    else if (nstep>0) {
      ierr = PVodeReInit(gdata);CHKERRQ(ierr);
    }
  }

  /* renormalize M _after_ WriteLog to have consistent data in logfile */

  PetscReal devNorm;
  devNorm=RenormVec(gdata->M,1.0,renormtol,ND);

  if (PetscAbsReal(devNorm) > renormtol) {
    PetscPrintf(PETSC_COMM_WORLD,
      "t=%g ns: renormalized M because |%g| > %g\n",
      gdata->time*gdata->tscale*1e9,
      devNorm,
      renormtol
    );

    /* if renormalization occurs within 3 timesteps */
    if (nstep<3 && reltol > 1e-10) {
      /* adjust tolerances */
      ierr = adjtol(gdata,1.0/tolscale);CHKERRQ(ierr);
    }
    else {
      /* only reinit PVODE after renormalization */
      ierr = PVodeReInit(gdata);CHKERRQ(ierr);
    }
  }
  /* check if timestep is very small */
  else if ((gdata->time-timelast)*gdata->tscale*1e9<1e-5 && nstep>100) {
    PetscInfo2(0,
      "Warning (t=%g ns): timestep is %g < 1e-4\n",
      gdata->time*gdata->tscale*1e9,
      (gdata->time-timelast)*gdata->tscale*1e9
    );
/*
    ierr = adjtol(gdata,1.0/tolscale);CHKERRQ(ierr);
    ierr = adjtoldem(gdata,1.0/tolscale);CHKERRQ(ierr);
*/

    if (gdata->VHdem!=PETSC_NULL &&
          (gdata->Edem-Edemold)*dEdem1<0 &&
          (gdata->Edem-Edemold)*dEdem2>0) {
      PetscInfo2(0,
        "and Edem oscillates: %g %g\n",
        gdata->Edem,
        Edemold
      );
/*
      PetscPrintf(PETSC_COMM_WORLD,
        "and Edem oscillates: %g %g %g %g\n",
        gdata->Edem,
        Edemold,
        gdata->Edem-Edemold,
        dEdem1,
        dEdem2
      );
*/
/*
      ierr = adjtoldem(gdata,1.0/tolscale);CHKERRQ(ierr);
*/
    }
  }
  /* check if Edem oscillates */
  else if (gdata->VHdem!=PETSC_NULL &&
           (gdata->Edem-Edemold)*dEdem1<0 &&
           (gdata->Edem-Edemold)*dEdem2>0) {
    PetscInfo2(0,
      "Warning (t=%g ns): Edem oscillates: %g\n",
      gdata->time*gdata->tscale*1e9,
      gdata->Edem
    );
/*
    PetscPrintf(PETSC_COMM_WORLD,
      "Warning (t=%g ns): Edem oscillates: %g %g %g %g %g\n",
      gdata->time*gdata->tscale*1e9,
      gdata->Edem,
      Edemold,
      gdata->Edem-Edemold,
      dEdem1,
      dEdem2
    );
*/
/*
    ierr = adjtoldem(gdata,1.0/tolscale);CHKERRQ(ierr);
*/
  }
  /* don't do this !?
    may lead to
    large tolerances->large error->norm not preserved->
    minimum time step->very slow simulation
  */
  else {
    long int ncfn;
    extern void *cvode_mem;
    CVodeGetNumNonlinSolvConvFails(cvode_mem,&ncfn);
    if (nstep > 30 &&
        devNorm <= renormtol/100.0 &&
        reltol <= tol/10.0 &&
        reltol <= 1e-2 &&
        ncfn == 0) {
      ierr = adjtol(gdata,tolscale);CHKERRQ(ierr);
    }
  }

  timelast=gdata->time;

  dEdem2=dEdem1;
  dEdem1=gdata->Edem-Edemold;
  Edemold=gdata->Edem;

  MagparFunctionProfReturn(0);
}