checkiterationllg.c

Go to the documentation of this file.
00001 /*
00002     This file is part of magpar.
00003 
00004     Copyright (C) 2002-2010 Werner Scholz
00005 
00006     www:   http://www.magpar.net/
00007     email: magpar(at)magpar.net
00008 
00009     magpar is free software; you can redistribute it and/or modify
00010     it under the terms of the GNU General Public License as published by
00011     the Free Software Foundation; either version 2 of the License, or
00012     (at your option) any later version.
00013 
00014     magpar is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00017     GNU General Public License for more details.
00018 
00019     You should have received a copy of the GNU General Public License
00020     along with magpar; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 */
00023 
00024 static char const Id[] = "$Id: checkiterationllg.c 2986 2010-02-23 18:20:01Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";
00026 
00027 #include "llg.h"
00028 #include "io/magpario.h"
00029 #include "util/util.h"
00030 #include "field/field.h"
00031 
00032 #ifdef ADDONS
00033 #include "addons/addons.h"
00034 #endif
00035 
00036 static int         doinit=1;
00037 
00038 static PetscReal   tol;
00039 static PetscReal   renormtol=-1.0;
00040 static int         condinp_equil;
00041 static PetscReal   condinp_j;
00042 static PetscReal   condinp_equil_j;
00043 static PetscReal   condinp_t;
00044 static PetscReal   *condinp_user=PETSC_NULL;
00045 static int         userinptot = 0;
00046 
00047 static PetscRandom rctx=PETSC_NULL;
00048 static PetscReal   magdist = 0.0;
00049 static PetscReal   magdist_ns = 1.0;
00050 static PetscReal   magdist_last_time = PETSC_MIN;
00051 
00052 /* TODO: resolve this extern dependency */
00053 extern KSP ksp_Au1;
00054 extern KSP ksp_Au2;
00055 extern realtype abstol, reltol; 
00057 int adjtoldem(GridData *gdata, double fac)
00058 {
00059   MagparFunctionInfoBegin;
00060 
00061   if (gdata->VHdem==PETSC_NULL) {
00062     MagparFunctionInfoReturn(0);
00063   }
00064 
00065   PetscReal rtol, atol, dtol;
00066   int maxits;
00067 
00068   ierr = KSPGetTolerances(ksp_Au1,&rtol,&atol,&dtol,(PetscInt*)&maxits);CHKERRQ(ierr);
00069 
00070   if (rtol < 1e-10) {
00071     MagparFunctionInfoReturn(0);
00072   }
00073 
00074   atol *= fac;
00075   rtol *= fac;
00076   ierr = KSPSetTolerances(ksp_Au1,rtol,atol,dtol,maxits);CHKERRQ(ierr);
00077 
00078   PetscReal time_ns = gdata->time*gdata->tscale*1e9;
00079   PetscPrintf(PETSC_COMM_WORLD,
00080     "at t=%g ns: ",time_ns
00081   );
00082 
00083   PetscPrintf(PETSC_COMM_WORLD,
00084     "%s KSP reltol to %g , abstol to %g\n",
00085     fac>1.0 ? "Increased" : "Reduced",
00086     rtol,atol
00087   );
00088 
00089   ierr = KSPGetTolerances(ksp_Au2,&rtol,&atol,&dtol,(PetscInt*)&maxits);CHKERRQ(ierr);
00090   atol *= fac;
00091   rtol *= fac;
00092   ierr = KSPSetTolerances(ksp_Au2,rtol,atol,dtol,maxits);CHKERRQ(ierr);
00093 
00094   MagparFunctionInfoReturn(0);
00095 }
00096 
00097 
00098 int adjtol(GridData *gdata, double fac)
00099 {
00100   MagparFunctionInfoBegin;
00101 
00102   if (reltol < 1e-10) {
00103     MagparFunctionInfoReturn(0);
00104   }
00105 
00106   PetscReal time_ns = gdata->time*gdata->tscale*1e9;
00107   PetscPrintf(PETSC_COMM_WORLD,
00108     "at t=%g ns: ",time_ns
00109   );
00110 
00111   adjtoldem(gdata,fac);
00112 
00113   reltol *= fac;
00114   abstol *= fac;
00115 
00116   extern void *cvode_mem;
00117   long int nstep;
00118   int flag;
00119   flag=CVodeGetNumSteps(cvode_mem,&nstep);
00120   PetscPrintf(PETSC_COMM_WORLD,
00121     "%s PVODE reltol to %g , abstol to %g ; nstep since last reinit: %i\n",
00122     fac>1.0 ? "Increased" : "Reduced",
00123     reltol,abstol,nstep
00124   );
00125 
00126   ierr = PVodeReInit(gdata);CHKERRQ(ierr);
00127 
00128   MagparFunctionInfoReturn(0);
00129 }
00130 
00131 
00132 int EquilCheck(Vec M,PetscReal time,PetscReal *torque)
00133 {
00134   static PetscReal tlast=time-0.00001;
00135   static Vec Mold=PETSC_NULL;
00136 
00137   MagparFunctionInfoBegin;
00138 
00139   /* initialize upon first call */
00140   if (Mold==PETSC_NULL) {
00141     ierr = VecDuplicate(M,&Mold);CHKERRQ(ierr);
00142     ierr = VecZeroEntries(Mold);CHKERRQ(ierr);
00143 
00144     /* copy M */
00145     ierr = VecCopy(M,Mold);CHKERRQ(ierr);
00146     /* return fake value since we have not done any timesteps yet
00147        just larger than the tolerance to keep the simulation running
00148     */
00149     *torque=tol*1.01;
00150 
00151     MagparFunctionInfoReturn(0);
00152   }
00153 
00154   ierr = VecAXPY(Mold,-1.0,M);CHKERRQ(ierr);
00155 
00156   /* calculate norm of delta M */
00157   *torque=RenormVec(Mold,0.0,PETSC_MAX,ND)/(time-tlast);
00158 
00159   ierr = VecCopy(M,Mold);CHKERRQ(ierr);
00160 
00161   tlast=time;
00162 
00163   MagparFunctionProfReturn(0);
00164 }
00165 
00166 
00167 int CheckIterationLLG_Init(GridData *gdata)
00168 {
00169   PetscTruth flg;
00170 
00171   MagparFunctionLogBegin;
00172   char fn[256];
00173   char line[256];
00174   int rank;
00175   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00176 
00177   ierr = PetscOptionsGetReal(PETSC_NULL,"-renormtol",&renormtol,&flg);CHKERRQ(ierr);
00178   if (flg!=PETSC_TRUE) {
00179     renormtol=1e-2;
00180     PetscPrintf(PETSC_COMM_WORLD,
00181       "Option -renormtol not found, using default value: %g\n",
00182       renormtol
00183     );
00184   }
00185   PetscPrintf(PETSC_COMM_WORLD,"renormtol:   %g\n",renormtol);
00186 
00187   ierr = PetscOptionsGetReal(PETSC_NULL,"-tol",&tol,&flg);CHKERRQ(ierr);
00188   if (flg!=PETSC_TRUE) {
00189     tol=1e-5;
00190     PetscPrintf(PETSC_COMM_WORLD,
00191       "Option -tol not found, using default value: %g\n",
00192       tol
00193     );
00194   }
00195   PetscPrintf(PETSC_COMM_WORLD,"tol (equilibrium): %g\n",tol);
00196 
00197   if (abstol > tol/10.0) {
00198     PetscPrintf(PETSC_COMM_WORLD,
00199       "Warning: abstol=%g should be < (tol=%g)/10.0!\n",
00200       abstol,
00201       tol
00202     );
00203   }
00204 
00205   ierr = PetscOptionsGetInt(PETSC_NULL,"-condinp_equil", (PetscInt*)&condinp_equil, &flg);CHKERRQ(ierr);
00206   if (flg!=PETSC_TRUE) {
00207     condinp_equil=1;
00208     PetscPrintf(PETSC_COMM_WORLD,
00209       "Option -condinp_equil not found, using default value: %i\n",
00210       condinp_equil
00211     );
00212   }
00213   ierr = PetscOptionsGetReal(PETSC_NULL,"-condinp_j",&condinp_j,&flg);CHKERRQ(ierr);
00214   if (flg!=PETSC_TRUE) {
00215     condinp_j=0.1;
00216     PetscPrintf(PETSC_COMM_WORLD,
00217       "Option -condinp_j not found, using default value: %g\n",
00218       condinp_j
00219     );
00220   }
00221   ierr = PetscOptionsGetReal(PETSC_NULL,"-condinp_equil_j",&condinp_equil_j,&flg);CHKERRQ(ierr);
00222   if (flg!=PETSC_TRUE) {
00223     condinp_equil_j=0.0;
00224     PetscPrintf(PETSC_COMM_WORLD,
00225       "Option -condinp_equil_j not found, using default value: %g\n",
00226       condinp_equil_j
00227     );
00228   }
00229   ierr = PetscOptionsGetReal(PETSC_NULL,"-condinp_t",&condinp_t,&flg);CHKERRQ(ierr);
00230   if (flg!=PETSC_TRUE) {
00231     condinp_t=PETSC_MAX;
00232     PetscPrintf(PETSC_COMM_WORLD,
00233       "Option -condinp_t not found, using default value: %g\n",
00234       condinp_t
00235     );
00236   }
00237 
00238   /* user times in file name specified by -condinp_user */
00239   ierr = PetscOptionsGetString(PETSC_NULL,"-condinp_file_t_ns",fn,sizeof(fn),&flg);CHKERRQ(ierr);
00240   if (flg!=PETSC_TRUE) {
00241     PetscPrintf(PETSC_COMM_WORLD,
00242       "Option -condinp_file_t_ns not found, user inp output times NOT activated\n"
00243     );
00244     userinptot = 0;
00245   }
00246   else {
00247     if (!rank) {
00248       FILE *fd=NULL;
00249       PetscPrintf(PETSC_COMM_WORLD,"Opening condinp_file_t_ns file %s\n",fn);
00250       ierr = PetscFOpen(PETSC_COMM_WORLD,fn,"r",&fd);CHKERRQ(ierr);
00251       if (!fd){
00252         SETERRQ1(PETSC_ERR_FILE_OPEN,"Cannot open %s\n",fn);
00253       }
00254       read_one_line(fd,line,sizeof(line),NULL);
00255       if (sscanf(line,"%i",&userinptot) != 1) {
00256         SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"could not read number of user inp values\n");
00257       }
00258       /* alloc memory to store user time values */
00259       ierr = PetscMalloc(userinptot*sizeof(PetscReal),&condinp_user);CHKERRQ(ierr);
00260 
00261       for (int i=0;i<userinptot;i++) {
00262         read_one_line(fd,line,sizeof(line),NULL);
00263         if (sscanf(line,"%lf",condinp_user + i) != 1) {
00264         SETERRQ1(PETSC_ERR_FILE_UNEXPECTED,"could not read user inp time value %i\n",i);
00265         }
00266         condinp_user[i] /= 1e9*gdata->tscale;
00267       }
00268 
00269       PetscPrintf(PETSC_COMM_WORLD,"Number of user inp output times read in: %i\n",userinptot);
00270       ierr = PetscFClose(PETSC_COMM_WORLD,fd);CHKERRQ(ierr);
00271     }
00272     /* broadcast condinp_user data to all processors */
00273     ierr = MPI_Bcast(&userinptot,1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
00274     if (rank) {
00275       ierr = PetscMalloc(userinptot*sizeof(PetscReal),&condinp_user);CHKERRQ(ierr);
00276     }
00277     ierr = MPI_Bcast(condinp_user,userinptot,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
00278   }
00279   condinp_t /= 1e9*gdata->tscale;
00280 
00281   /* disturb magnetization distribution with random vector */
00282   ierr = PetscOptionsGetReal(PETSC_NULL,"-magdist_llg", &magdist,&flg);CHKERRQ(ierr);
00283   if (flg!=PETSC_TRUE) {
00284     magdist=0.0;
00285     PetscPrintf(PETSC_COMM_WORLD,
00286       "Option -magdist_llg not found, using default value: %g\n",
00287       magdist
00288     );
00289   }
00290 
00291   if (magdist>0.0) {
00292     ierr = PetscOptionsGetReal(PETSC_NULL,"-magdist_ns", &magdist_ns,&flg);CHKERRQ(ierr);
00293     if (flg!=PETSC_TRUE) {
00294       magdist_ns=1.0;
00295       PetscPrintf(PETSC_COMM_WORLD,
00296         "Option -magdist_ns not found, using default value: %g\n",
00297         magdist_ns
00298       );
00299     }
00300 
00301     PetscReal time_ns = gdata->time*gdata->tscale*1e9;
00302     magdist_last_time = time_ns;
00303     PetscPrintf(PETSC_COMM_WORLD,
00304       "Initializing magdist_last_time: %g ns\n",
00305       magdist_last_time
00306     );
00307 
00308     ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
00309     ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
00310 
00311     PetscPrintf(PETSC_COMM_WORLD,"Applying magnetization disturbance of %g every %g ns.\n",magdist,magdist_ns);
00312   }
00313 
00314   MagparFunctionLogReturn(0);
00315 }
00316 
00317 
00318 int CheckIterationLLG(GridData *gdata)
00319 {
00320   static PetscReal Etotold=PETSC_MAX;
00321   static PetscReal hextold=PETSC_MAX;
00322   static PetscReal condinp_last_t=PETSC_MIN;
00323   static PetscReal condinp_last_j=PETSC_MIN;
00324   static PetscReal condinp_last_equil_j=PETSC_MIN;
00325 
00326   const PetscReal  tolscale=3.0;
00327 
00328   static int nsteps=-1;
00329   static int csteps;
00330 
00331   static PetscReal Edemold=gdata->Edem;
00332   static PetscReal dEdem1=0;
00333   static PetscReal dEdem2=0;
00334 
00335   static PetscReal timelast=gdata->time;
00336 
00337   static int userinpindex = 0;
00338 
00339 #undef CHKDMDT
00340 #ifdef CHKDMDT
00341   static Vec dMdt;
00342 #endif
00343 
00344   MagparFunctionInfoBegin;
00345 
00346   PetscReal time_ns = gdata->time*gdata->tscale*1e9;
00347 
00348   /* initialize if necessary */
00349   if (doinit) {
00350     ierr = CheckIterationLLG_Init(gdata);CHKERRQ(ierr);
00351 #ifdef CHKDMDT
00352     VecDuplicate(gdata->M,&dMdt);
00353 #endif
00354     doinit=0;
00355   }
00356 
00357 #ifdef CHKDMDT
00358   /* This method could be more expensive than the old method because
00359      we recalculate all the fields when we call calc_dMdt!
00360      However, we would get access the old dMdt - then it would be cheaper.
00361   */
00362   PetscReal veq;
00363   calc_dMdt(PETSC_NULL,0.0,gdata->M,dMdt,gdata);
00364   veq=RenormVec(dMdt,0.0,PETSC_MAX,ND);
00365   gdata->equil=veq;
00366 
00367 #if 1
00368   ierr = EquilCheck(gdata->M,gdata->time,&gdata->vequil);CHKERRQ(ierr);
00369 #endif
00370 
00371   PetscPrintf(PETSC_COMM_WORLD,"torque = %g  %g  %g\n",veq,gdata->vequil,veq/gdata->vequil);
00372 #else
00373   ierr = EquilCheck(gdata->M,gdata->time,&gdata->vequil);CHKERRQ(ierr);
00374 #endif
00375 
00376   /* We could also use differential dMdt calculated in calc_dMdt through RHSfunction (see there).
00377      However, there is more noise than if we take the finite difference over one full time step.
00378   */
00379 
00380   if (gdata->vequil < tol) {
00381     /* count how long we have been in equilibrium in succession */
00382     gdata->equil++;
00383   }
00384   else {
00385     /* reset counter if we are not in equilibrium */
00386     gdata->equil=0;
00387   }
00388 
00389   /* check if total energy increased (which it should not)
00390      print warning only if external field has not changed
00391   */
00392   PetscReal hext;
00393   hext=Hexternal_hext();
00394   if (gdata->Etot > Etotold &&
00395       hext == hextold &&
00396       ((gdata->Etot-Etotold)/Etotold)>1e-8) {
00397 
00398     PetscPrintf(MPI_COMM_WORLD,
00399       "Warning (t=%g ns): Etot increased by %g (%g*Etot) from %g to %g J/m^3\n",
00400       time_ns,
00401       (gdata->Etot-Etotold)*gdata->escale,
00402       (gdata->Etot-Etotold)/Etotold,
00403       Etotold*gdata->escale,
00404       gdata->Etot*gdata->escale
00405     );
00406 
00407 /*
00408     ierr = adjtol(gdata,1.0/tolscale);CHKERRQ(ierr);
00409 */
00410   }
00411   Etotold=gdata->Etot;
00412   hextold=hext;
00413 
00414   PetscReal valMpHext;
00415   ierr = MpHext(gdata->M,&valMpHext);CHKERRQ(ierr);
00416 
00417   /* check if a set of output data should be written */
00418   if (PetscAbsReal(valMpHext-condinp_last_j) > condinp_j) {
00419     gdata->inp=PetscAbsInt(gdata->inp);
00420 
00421     condinp_last_t=gdata->time;
00422     condinp_last_j=valMpHext;
00423   }
00424   else if (PetscAbsReal(gdata->time-condinp_last_t) > condinp_t) {
00425     gdata->inp=PetscAbsInt(gdata->inp);
00426 
00427     condinp_last_t=gdata->time;
00428     condinp_last_j=valMpHext;
00429   }
00430   else if ( userinpindex < userinptot && gdata->time >= condinp_user[userinpindex]) {
00431       gdata->inp=PetscAbsInt(gdata->inp);
00432 
00433       condinp_last_t=gdata->time;
00434       condinp_last_j=valMpHext;
00435       userinpindex++;
00436   }
00437   else if (condinp_equil==1 && gdata->equil > 0) {
00438     /* check for condinp_equil_j (defaults to 0.0) */
00439     if (PetscAbsReal(valMpHext-condinp_last_equil_j) > condinp_equil_j) {
00440       gdata->inp=PetscAbsInt(gdata->inp);
00441 
00442       condinp_last_t=gdata->time;
00443       condinp_last_j=valMpHext;
00444       condinp_last_equil_j=valMpHext;
00445     }
00446   }
00447 
00448   if (nsteps<0) {
00449     PetscTruth flg;
00450     ierr = PetscOptionsGetInt(PETSC_NULL,"-ts_logsteps",(PetscInt*)&nsteps,&flg);CHKERRQ(ierr);
00451     if (flg!=PETSC_TRUE) {
00452       nsteps=1;
00453       PetscPrintf(PETSC_COMM_WORLD,
00454         "Option -ts_logsteps not found, using default value: %i\n",
00455         nsteps
00456       );
00457     }
00458     if (nsteps<1) {
00459       SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Option set: -ts_logsteps %i, but has to be >0\n",nsteps);
00460     }
00461     csteps=0;
00462   }
00463 
00464   if (csteps % nsteps == 0 || gdata->equil>0) {
00465     ierr = WriteLog(gdata);CHKERRQ(ierr);
00466     ierr = WriteSet(gdata);CHKERRQ(ierr);
00467     csteps=0;
00468   }
00469   csteps++;
00470 
00471 #ifdef ADDONS
00472   /* move magnetization if requested */
00473   ierr = MovingField(gdata);CHKERRQ(ierr);
00474 #endif
00475 
00476   extern void *cvode_mem;
00477   long int nstep;
00478   int flag;
00479   flag=CVodeGetNumSteps(cvode_mem,&nstep);
00480 
00481   /* update external field if system is in equilibrium */
00482   if (gdata->equil>0) {
00483     ierr = Hext_ho_hstep(gdata);CHKERRQ(ierr);
00484 
00485     /* check if we are still in equilibrium */
00486     if (gdata->equil>0) {
00487       /* that's because the field has not been changed */
00488       PetscPrintf(PETSC_COMM_WORLD,"gdata->equil %i >= 1\n",gdata->equil);
00489       gdata->mode=99;
00490     }
00491     /* restart PVODE because the external field is changed abruptly */
00492     else if (nstep>0) {
00493       ierr = PVodeReInit(gdata);CHKERRQ(ierr);
00494     }
00495   }
00496 
00497 
00498   /* disturb magnetization distribution with random vector */
00499   if (magdist>0.0) {
00500     if (magdist_last_time==PETSC_MIN) {
00501       magdist_last_time = time_ns;
00502     }
00503     else if (time_ns - magdist_last_time > magdist_ns) {
00504       /* randomly distort M */
00505       PetscPrintf(PETSC_COMM_WORLD,"Disturbing magnetization distribution at %g ns by %g.\n",time_ns,magdist);
00506       ierr = DistortVec(rctx,gdata->M,magdist);CHKERRQ(ierr);
00507 
00508       /* renormalize M */
00509       PetscReal devNorm;
00510       devNorm=RenormVec(gdata->M,1.0,D_EPS,ND);
00511 
00512       ierr = PVodeReInit(gdata);CHKERRQ(ierr);
00513 
00514       magdist_last_time=time_ns;
00515     }
00516   }
00517 
00518 
00519   /* renormalize M _after_ WriteLog to have consistent data in logfile */
00520 
00521   PetscReal devNorm;
00522   devNorm=RenormVec(gdata->M,1.0,renormtol,ND);
00523 
00524   if (PetscAbsReal(devNorm) > renormtol) {
00525     PetscPrintf(PETSC_COMM_WORLD,
00526       "t=%g ns: renormalized M because |%g| > %g\n",
00527       time_ns,
00528       devNorm,
00529       renormtol
00530     );
00531 
00532     /* if renormalization occurs within 3 timesteps */
00533     if (nstep<3 && reltol > 1e-10) {
00534       /* adjust tolerances */
00535       ierr = adjtol(gdata,1.0/tolscale);CHKERRQ(ierr);
00536     }
00537     else {
00538       /* only reinit PVODE after renormalization */
00539       ierr = PVodeReInit(gdata);CHKERRQ(ierr);
00540     }
00541   }
00542   /* check if timestep is very small */
00543   else if ((gdata->time-timelast)*gdata->tscale*1e9<1e-5 && nstep>100) {
00544     PetscInfo2(0,
00545       "Warning (t=%g ns): timestep is %g < 1e-4\n",
00546       time_ns,
00547       (gdata->time-timelast)*gdata->tscale*1e9
00548     );
00549 /*
00550     ierr = adjtol(gdata,1.0/tolscale);CHKERRQ(ierr);
00551     ierr = adjtoldem(gdata,1.0/tolscale);CHKERRQ(ierr);
00552 */
00553 
00554     if (gdata->VHdem!=PETSC_NULL &&
00555           (gdata->Edem-Edemold)*dEdem1<0 &&
00556           (gdata->Edem-Edemold)*dEdem2>0) {
00557       PetscInfo2(0,
00558         "and Edem oscillates: %g %g\n",
00559         gdata->Edem,
00560         Edemold
00561       );
00562 /*
00563       PetscPrintf(PETSC_COMM_WORLD,
00564         "and Edem oscillates: %g %g %g %g\n",
00565         gdata->Edem,
00566         Edemold,
00567         gdata->Edem-Edemold,
00568         dEdem1,
00569         dEdem2
00570       );
00571 */
00572 /*
00573       ierr = adjtoldem(gdata,1.0/tolscale);CHKERRQ(ierr);
00574 */
00575     }
00576   }
00577   /* check if Edem oscillates */
00578   else if (gdata->VHdem!=PETSC_NULL &&
00579            (gdata->Edem-Edemold)*dEdem1<0 &&
00580            (gdata->Edem-Edemold)*dEdem2>0) {
00581     PetscInfo2(0,
00582       "Warning (t=%g ns): Edem oscillates: %g\n",
00583       time_ns,
00584       gdata->Edem
00585     );
00586 /*
00587     PetscPrintf(PETSC_COMM_WORLD,
00588       "Warning (t=%g ns): Edem oscillates: %g %g %g %g %g\n",
00589       time_ns,
00590       gdata->Edem,
00591       Edemold,
00592       gdata->Edem-Edemold,
00593       dEdem1,
00594       dEdem2
00595     );
00596 */
00597 /*
00598     ierr = adjtoldem(gdata,1.0/tolscale);CHKERRQ(ierr);
00599 */
00600   }
00601   /* don't do this !?
00602     may lead to
00603     large tolerances->large error->norm not preserved->
00604     minimum time step->very slow simulation
00605   */
00606   else {
00607     long int ncfn;
00608     extern void *cvode_mem;
00609     CVodeGetNumNonlinSolvConvFails(cvode_mem,&ncfn);
00610     if (nstep > 30 &&
00611         devNorm <= renormtol/100.0 &&
00612         reltol <= tol/10.0 &&
00613         reltol <= 1e-2 &&
00614         ncfn == 0) {
00615       ierr = adjtol(gdata,tolscale);CHKERRQ(ierr);
00616     }
00617   }
00618 
00619   timelast=gdata->time;
00620 
00621   dEdem2=dEdem1;
00622   dEdem1=gdata->Edem-Edemold;
00623   Edemold=gdata->Edem;
00624 
00625   MagparFunctionProfReturn(0);
00626 }

magpar - Parallel Finite Element Micromagnetics Package
Copyright (C) 2002-2010 Werner Scholz