00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
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
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
00140 if (Mold==PETSC_NULL) {
00141 ierr = VecDuplicate(M,&Mold);CHKERRQ(ierr);
00142 ierr = VecZeroEntries(Mold);CHKERRQ(ierr);
00143
00144
00145 ierr = VecCopy(M,Mold);CHKERRQ(ierr);
00146
00147
00148
00149 *torque=tol*1.01;
00150
00151 MagparFunctionInfoReturn(0);
00152 }
00153
00154 ierr = VecAXPY(Mold,-1.0,M);CHKERRQ(ierr);
00155
00156
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
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
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
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
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
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
00359
00360
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
00377
00378
00379
00380 if (gdata->vequil < tol) {
00381
00382 gdata->equil++;
00383 }
00384 else {
00385
00386 gdata->equil=0;
00387 }
00388
00389
00390
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
00409
00410 }
00411 Etotold=gdata->Etot;
00412 hextold=hext;
00413
00414 PetscReal valMpHext;
00415 ierr = MpHext(gdata->M,&valMpHext);CHKERRQ(ierr);
00416
00417
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
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
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
00482 if (gdata->equil>0) {
00483 ierr = Hext_ho_hstep(gdata);CHKERRQ(ierr);
00484
00485
00486 if (gdata->equil>0) {
00487
00488 PetscPrintf(PETSC_COMM_WORLD,"gdata->equil %i >= 1\n",gdata->equil);
00489 gdata->mode=99;
00490 }
00491
00492 else if (nstep>0) {
00493 ierr = PVodeReInit(gdata);CHKERRQ(ierr);
00494 }
00495 }
00496
00497
00498
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
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
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
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
00533 if (nstep<3 && reltol > 1e-10) {
00534
00535 ierr = adjtol(gdata,1.0/tolscale);CHKERRQ(ierr);
00536 }
00537 else {
00538
00539 ierr = PVodeReInit(gdata);CHKERRQ(ierr);
00540 }
00541 }
00542
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
00551
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
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575 }
00576 }
00577
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
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600 }
00601
00602
00603
00604
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 }