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: eminisolve.c 3012 2010-03-26 16:05:11Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027 #include "emini.h"
00028 #include "field/field.h"
00029 #include "util/util.h"
00030 #include "llg/llg.h"
00031 #include "io/magpario.h"
00032
00033 #ifdef ADDONS
00034 #include "addons/addons.h"
00035 #endif
00036
00037 #undef DEBUG
00038 #ifdef DEBUG
00039 static Vec dMdt=PETSC_NULL;
00040 #endif
00041
00042
00043
00044
00045 #undef DAMP
00046
00047
00048
00049
00050
00051
00052 static TAO_SOLVER tao;
00053 static TAO_APPLICATION taoapp;
00054 static PetscRandom rctx=PETSC_NULL;
00055 static PetscReal magdist;
00056 static PetscReal tol;
00057 static int tao_max_its;
00058
00059 static Vec Msphere;
00060 static Vec VHtotsphere;
00061
00062 #if 0
00063 #define NO_VMS3
00064 static Vec X2;
00065 #endif
00066
00067 #undef TESTTORQ
00068 #ifdef TESTTORQ
00069 static Vec dMdt=PETSC_NULL;
00070 #endif
00071
00072 #ifdef DAMP
00073 static Vec grad;
00074 #endif
00075
00076 int TaoEvalEnergyGrad(TAO_APPLICATION taoapp,Vec X,double *f, Vec G,void *ptr)
00077 {
00078 GridData *gdata = (GridData*)ptr;
00079
00080 MagparFunctionInfoBegin;
00081
00082 if (gdata->mode==53) {
00083
00084
00085
00086
00087
00088 ierr = VecCopy(X,gdata->M);CHKERRQ(ierr);
00089 }
00090 else {
00091 #ifdef NO_VMS3
00092 ierr = VecCopy(X,X2);CHKERRQ(ierr);
00093
00094 PetscReal *ta_X, *ta_VMs3;
00095 ierr = VecGetArray(X2,&ta_X);CHKERRQ(ierr);
00096 ierr = VecGetArray(gdata->VMs3,&ta_VMs3);CHKERRQ(ierr);
00097 for (int i=0; i<gdata->ln_vert; i++) {
00098 ta_X[2*i+0] /= ta_VMs3[ND*i];
00099 ta_X[2*i+1] /= ta_VMs3[ND*i];
00100 }
00101 ierr = VecRestoreArray(X2,&ta_X);CHKERRQ(ierr);
00102 ierr = VecRestoreArray(gdata->VMs3,&ta_VMs3);CHKERRQ(ierr);
00103
00104 ierr = Sphere2Cart(X2,gdata->M);CHKERRQ(ierr);
00105 #else
00106 ierr = Sphere2Cart(X,gdata->M);CHKERRQ(ierr);
00107 #endif
00108 }
00109
00110 ierr = Htot_Gradient(gdata);CHKERRQ(ierr);
00111
00112 #ifdef DAMP
00113
00114
00115 PetscReal *ta_M, *ta_Htot, *ta_grad;
00116 ierr = VecGetArray(gdata->M,&ta_M);CHKERRQ(ierr);
00117 ierr = VecGetArray(gdata->VHtot,&ta_Htot);CHKERRQ(ierr);
00118 ierr = VecGetArray(grad,&ta_grad);CHKERRQ(ierr);
00119
00120 for (int i=0; i<gdata->ln_vert; i++) {
00121 PetscReal *damp, mxh[ND];
00122 PetscReal a1;
00123
00124
00125 if (gdata->propdat[NP*gdata->vertprop[i]+4]<=0.0) continue;
00126
00127 damp=ta_grad+ND*i;
00128
00129 douter(ND,ta_M+ND*i,ta_Htot+ND*i,mxh);
00130 douter(ND,ta_M+ND*i,mxh,damp);
00131
00132
00133 a1=0.001;
00134 my_dscal(ND,a1,damp,1);
00135 }
00136
00137 ierr = VecRestoreArray(gdata->M,&ta_M);CHKERRQ(ierr);
00138 ierr = VecRestoreArray(gdata->VHtot,&ta_Htot);CHKERRQ(ierr);
00139 ierr = VecRestoreArray(grad,&ta_grad);CHKERRQ(ierr);
00140 #endif
00141
00142 ierr = Htot_Energy(gdata);CHKERRQ(ierr);
00143 *f=gdata->Etot;
00144
00145
00146
00147
00148
00149
00150
00151 #ifdef DAMP
00152 if (gdata->mode==53) {
00153 ierr = VecCopy(grad,G);CHKERRQ(ierr);
00154 }
00155 else {
00156 ierr = Cart2SphereDiff(grad,X,G);CHKERRQ(ierr);
00157 }
00158 #else
00159 if (gdata->mode==53) {
00160 ierr = VecCopy(gdata->VHtot,G);CHKERRQ(ierr);
00161 }
00162 else {
00163 ierr = Cart2SphereDiff(gdata->VHtot,X,G);CHKERRQ(ierr);
00164
00165 #ifdef NO_VMS3
00166 PetscReal *ta_X, *ta_VMs3;
00167 ierr = VecGetArray(G,&ta_X);CHKERRQ(ierr);
00168 ierr = VecGetArray(gdata->VMs3,&ta_VMs3);CHKERRQ(ierr);
00169 for (int i=0; i<gdata->ln_vert; i++) {
00170 ta_X[2*i+0] *= ta_VMs3[ND*i];
00171 ta_X[2*i+1] *= ta_VMs3[ND*i];
00172 }
00173 ierr = VecRestoreArray(G,&ta_X);CHKERRQ(ierr);
00174 ierr = VecRestoreArray(gdata->VMs3,&ta_VMs3);CHKERRQ(ierr);
00175 #endif
00176 }
00177 #endif
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188 MagparFunctionProfReturn(0);
00189 }
00190
00191 int ConvTest(TAO_SOLVER tao,void *cctx)
00192 {
00193 GridData *gdata = (GridData*)cctx;
00194
00195 MagparFunctionInfoBegin;
00196
00197 int iter;
00198 double ff,gnorm,cnorm,xdiff;
00199 TaoTerminateReason reason;
00200 ierr = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,&cnorm,&xdiff,&reason);CHKERRQ(ierr);
00201 if (iter>tao_max_its) {
00202 ierr = TaoSetTerminationReason(tao,TAO_DIVERGED_MAXITS);CHKERRQ(ierr);
00203 MagparFunctionProfReturn(0);
00204 }
00205
00206 #ifdef TESTTORQ
00207
00208
00209
00210
00211 if (gdata->mode!=53) {
00212 Vec sol;
00213 ierr = TaoAppGetSolutionVec(taoapp,&sol);CHKERRQ(ierr);
00214 assert(sol==Msphere);
00215 ierr = Sphere2Cart(sol,gdata->M);CHKERRQ(ierr);
00216 }
00217
00218 ierr = calc_dMdt(PETSC_NULL,0.0,gdata->M,dMdt,gdata);CHKERRQ(ierr);
00219 gdata->vequil=RenormVec(dMdt,0.0,PETSC_MAX,ND);
00220
00221 ierr = TaoGetTerminationReason(tao,&reason);CHKERRQ(ierr);
00222
00223 PetscPrintf(PETSC_COMM_WORLD,"torque = %g Etot = %g %g reason = %i\n",
00224 gdata->vequil,gdata->Etot,gdata->Etot*gdata->escale,reason
00225 );
00226
00227 static int equil=0;
00228 if (gdata->vequil <= tol) {
00229 ierr = TaoSetTerminationReason(tao,TAO_CONVERGED_USER);CHKERRQ(ierr);
00230
00231 equil++;
00232 }
00233 else {
00234 ierr = TaoSetTerminationReason(tao,TAO_CONTINUE_ITERATING);CHKERRQ(ierr);
00235
00236 equil=0;
00237 }
00238
00239 #else
00240
00241 static PetscReal Etotold;
00242 static int equil;
00243
00244 if (iter==0) {
00245 equil=0;
00246 Etotold=PETSC_MAX;
00247 }
00248
00249 PetscLogDouble t_now;
00250 PetscGetTime(&t_now);
00251 static PetscLogDouble t_last=t_now;
00252
00253 PetscPrintf(PETSC_COMM_WORLD,"%5i %4i %6.2f %22.16e %22.16e %22.16e %22.16e %22.16e",
00254 iter,
00255 equil,
00256 t_now-t_last,
00257 gdata->Etot,gdata->Etot*gdata->escale,
00258 Etotold,Etotold*gdata->escale,
00259 (Etotold-gdata->Etot)/PetscAbsReal(gdata->Etot)
00260 );
00261 #ifdef DEBUG
00262 ierr = calc_dMdt(PETSC_NULL,0.0,gdata->M,dMdt,gdata);CHKERRQ(ierr);
00263 gdata->vequil=RenormVec(dMdt,0.0,PETSC_MAX,ND);
00264
00265 PetscPrintf(PETSC_COMM_WORLD,"%22.16e",
00266 gdata->vequil
00267 );
00268 #endif
00269 PetscPrintf(PETSC_COMM_WORLD,"\n");
00270 t_last=t_now;
00271
00272
00273
00274
00275
00276 #define NEQUIL 5
00277 if (PetscAbs(Etotold-gdata->Etot)/PetscAbsReal(gdata->Etot)<tol) {
00278 if (equil>=NEQUIL) {
00279 ierr = TaoSetTerminationReason(tao,TAO_CONVERGED_USER);CHKERRQ(ierr);
00280 equil=0;
00281 }
00282 else {
00283 ierr = TaoSetTerminationReason(tao,TAO_CONTINUE_ITERATING);CHKERRQ(ierr);
00284
00285 equil++;
00286 }
00287 }
00288 else {
00289 ierr = TaoSetTerminationReason(tao,TAO_CONTINUE_ITERATING);CHKERRQ(ierr);
00290
00291 equil=0;
00292 }
00293 Etotold=gdata->Etot;
00294 #endif
00295
00296 #ifdef DEBUG
00297 ierr = WriteLog(gdata);CHKERRQ(ierr);
00298 gdata->inp=PetscAbsInt(gdata->inp);
00299 ierr = WriteSet(gdata);CHKERRQ(ierr);
00300 #endif
00301
00302 MagparFunctionProfReturn(0);
00303 }
00304
00305
00306 int myTSCreateEmini(GridData *gdata)
00307 {
00308 MagparFunctionLogBegin;
00309
00310 gdata->time = 0.0;
00311
00312 #ifdef DAMP
00313 PetscPrintf(PETSC_COMM_WORLD,"DAMP defined: Using damping term for gradient/line search\n");
00314 ierr = VecDuplicate(gdata->M,&grad);CHKERRQ(ierr);
00315 #else
00316 PetscPrintf(PETSC_COMM_WORLD,"DAMP undefined: Using field for gradient/line search\n");
00317 #endif
00318
00319
00320 PetscTruth flg;
00321 int tao_lmm_vectors;
00322 ierr = PetscOptionsGetInt(PETSC_NULL,"-tao_lmm_vectors",(PetscInt*)&tao_lmm_vectors,&flg);CHKERRQ(ierr);
00323 if (flg!=PETSC_TRUE) {
00324 ierr = PetscOptionsSetValue("-tao_lmm_vectors","100");CHKERRQ(ierr);
00325 }
00326
00327 char str[256];
00328 PetscOptionsGetString(PETSC_NULL,"-tao_method",str,256,&flg);
00329 if (flg!=PETSC_TRUE) {
00330 ierr = PetscSNPrintf(str,255,"tao_lmvm");CHKERRQ(ierr);
00331 PetscPrintf(PETSC_COMM_WORLD,
00332 "Option -tao_method not found, using default value: %s\n",
00333 str
00334 );
00335 }
00336
00337 ierr = TaoCreate(PETSC_COMM_WORLD,str,&tao);CHKERRQ(ierr);
00338 ierr = TaoApplicationCreate(PETSC_COMM_WORLD,&taoapp);CHKERRQ(ierr);
00339 if (gdata->mode==53) {
00340 PetscPrintf(PETSC_COMM_WORLD,"mode %i: Using 3D Cartesian coordinates\n",gdata->mode);
00341
00342 PetscReal *ta_M;
00343 ierr = VecGetArray(gdata->M,&ta_M);CHKERRQ(ierr);
00344 for (int j=0;j<gdata->ln_vert;j++) {
00345 PetscReal ms;
00346 ms=gdata->propdat[NP*gdata->vertprop[j]+4];
00347
00348 if (ms==0.0) {
00349
00350
00351 ta_M[ND*j+0] = 0.0;
00352 ta_M[ND*j+1] = 0.0;
00353 ta_M[ND*j+2] = 0.0;
00354 }
00355 }
00356 ierr = VecRestoreArray(gdata->M,&ta_M);CHKERRQ(ierr);
00357
00358 ierr = TaoAppSetInitialSolutionVec(taoapp,gdata->M);CHKERRQ(ierr);
00359 }
00360 else {
00361 PetscPrintf(PETSC_COMM_WORLD,"Using spherical coordinates\n");
00362 ierr = VecCreate(PETSC_COMM_WORLD,&Msphere);CHKERRQ(ierr);
00363 ierr = VecSetSizes(Msphere,(ND-1)*gdata->ln_vert,(ND-1)*gdata->n_vert);CHKERRQ(ierr);
00364 ierr = VecSetFromOptions(Msphere);CHKERRQ(ierr);
00365
00366 ierr = VecDuplicate(Msphere,&VHtotsphere);CHKERRQ(ierr);
00367 #ifdef NO_VMS3
00368 ierr = VecDuplicate(Msphere,&X2);CHKERRQ(ierr);
00369 #endif
00370
00371 ierr = TaoAppSetInitialSolutionVec(taoapp,Msphere);CHKERRQ(ierr);
00372 }
00373 ierr = TaoAppSetObjectiveAndGradientRoutine(taoapp,TaoEvalEnergyGrad,(void *)gdata);CHKERRQ(ierr);
00374 ierr = TaoSetConvergenceTest(tao,ConvTest,(void *)gdata);CHKERRQ(ierr);
00375
00376
00377
00378
00379
00380
00381
00382 ierr = PetscOptionsHasName(PETSC_NULL,"-tao_ls_ftol",&flg);CHKERRQ(ierr);
00383 if (flg!=PETSC_TRUE) {
00384 ierr = PetscOptionsSetValue("-tao_ls_ftol","1e-20");CHKERRQ(ierr);
00385 }
00386 ierr = PetscOptionsHasName(PETSC_NULL,"-tao_ls_rtol",&flg);CHKERRQ(ierr);
00387 if (flg!=PETSC_TRUE) {
00388 ierr = PetscOptionsSetValue("-tao_ls_rtol","1e-1");CHKERRQ(ierr);
00389 }
00390 ierr = PetscOptionsHasName(PETSC_NULL,"-tao_ls_gtol",&flg);CHKERRQ(ierr);
00391 if (flg!=PETSC_TRUE) {
00392 ierr = PetscOptionsSetValue("-tao_ls_gtol","0.99");CHKERRQ(ierr);
00393 }
00394
00395 ierr = PetscOptionsGetInt(PETSC_NULL,"-tao_max_its",(PetscInt*)&tao_max_its,&flg);CHKERRQ(ierr);
00396 if (flg!=PETSC_TRUE) {
00397 tao_max_its=500;
00398 PetscPrintf(PETSC_COMM_WORLD,
00399 "Option -tao_max_its not found, using default value: %i\n",
00400 tao_max_its
00401 );
00402 }
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412 ierr = TaoSetOptions(taoapp,tao);CHKERRQ(ierr);
00413
00414
00415
00416
00417
00418
00419
00420 ierr = PetscOptionsGetReal(PETSC_NULL,"-magdist", &magdist,&flg);CHKERRQ(ierr);
00421 if (flg!=PETSC_TRUE) {
00422 magdist=0.0;
00423 PetscPrintf(PETSC_COMM_WORLD,
00424 "Option -magdist not found, using default value: %g\n",
00425 magdist
00426 );
00427 }
00428
00429 if (magdist<0.0) {
00430 PetscPrintf(PETSC_COMM_WORLD,"magdist = %g < 0.0: Distorting M by tilting in direction MxMxH\n",magdist);
00431 }
00432 else if (magdist>0.0) {
00433 PetscPrintf(PETSC_COMM_WORLD,"magdist = %g > 0.0: Distorting M randomly\n",magdist);
00434 }
00435 else {
00436 PetscPrintf(PETSC_COMM_WORLD,"magdist==0.0: Not distorting M\n");
00437 }
00438
00439 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
00440 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
00441
00442
00443 ierr = TaoView(tao);CHKERRQ(ierr);
00444
00445 ierr = PetscOptionsGetReal(PETSC_NULL,"-tol",&tol,&flg);CHKERRQ(ierr);
00446 if (flg!=PETSC_TRUE) {
00447 tol=1e-7;
00448 PetscPrintf(PETSC_COMM_WORLD,
00449 "Option -tol not found, using default value: %g\n",
00450 tol
00451 );
00452 }
00453 PetscPrintf(PETSC_COMM_WORLD,"tol (equilibrium): %g\n",tol);
00454
00455 #ifdef TESTTORQ
00456 ierr = VecDuplicate(gdata->M,&dMdt);
00457 #endif
00458 #ifdef DEBUG
00459 ierr = VecDuplicate(gdata->M,&dMdt);
00460 #endif
00461
00462 MagparFunctionLogReturn(0);
00463 }
00464
00465
00466 int EminiSolve(GridData *gdata)
00467 {
00468 MagparFunctionInfoBegin;
00469
00470 if (magdist<0.0) {
00471
00472 PetscReal *ta_M, *ta_Htot;
00473 VecGetArray(gdata->M,&ta_M);
00474 VecGetArray(gdata->VHtot,&ta_Htot);
00475
00476 for (int i=0; i<gdata->ln_vert; i++) {
00477 PetscReal damp[ND], mxh[ND];
00478
00479 douter(ND,ta_M+ND*i,ta_Htot+ND*i,mxh);
00480 douter(ND,ta_M+ND*i,mxh,damp);
00481
00482
00483 my_daxpy(ND,magdist,damp,1,ta_M+ND*i,1);
00484 }
00485
00486 VecRestoreArray(gdata->M,&ta_M);
00487 VecRestoreArray(gdata->VHtot,&ta_Htot);
00488 }
00489 else if (magdist>0.0) {
00490
00491 ierr = DistortVec(rctx,gdata->M,magdist);CHKERRQ(ierr);
00492
00493
00494 if (gdata->mode!=53) {
00495 PetscReal devNorm;
00496 devNorm=RenormVec(gdata->M,1.0,D_EPS,ND);
00497 }
00498 }
00499
00500 PetscPrintf(MPI_COMM_WORLD,"Starting energy minimization with TAO...\n");
00501 PetscPrintf(PETSC_COMM_WORLD,"iter equil CPUtstep (s) Etot (1) (J/m^3) Etotold (1) (J/m^3) tol torq\n");
00502 PetscPrintf(PETSC_COMM_WORLD,"--------------------------------------------------------------------------------------------------------------------------------------------------------------\n");
00503
00504 if (gdata->mode!=53) {
00505 ierr = Cart2Sphere(gdata->M,Msphere);CHKERRQ(ierr);
00506
00507 #ifdef NO_VMS3
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519 #endif
00520 }
00521
00522 ierr = TaoSolveApplication(taoapp,tao);CHKERRQ(ierr);
00523
00524 #ifdef ADDONS
00525
00526 #ifdef DEBUG
00527
00528 ierr = Htot(gdata);CHKERRQ(ierr);
00529 ierr = Htot_Energy(gdata);CHKERRQ(ierr);
00530
00531 gdata->inp=PetscAbsInt(gdata->inp);
00532 ierr = WriteLog(gdata);CHKERRQ(ierr);
00533 ierr = WriteSet(gdata);CHKERRQ(ierr);
00534 #endif
00535
00536 PetscTruth flg;
00537 int NSMOOTHSTEPS;
00538 ierr = PetscOptionsGetInt(PETSC_NULL,"-emini_mag_smooth_steps",(PetscInt*)&NSMOOTHSTEPS,&flg);CHKERRQ(ierr);
00539 if (flg!=PETSC_TRUE) {
00540 NSMOOTHSTEPS=1;
00541 }
00542
00543 if (NSMOOTHSTEPS>0) {
00544 PetscPrintf(MPI_COMM_WORLD,"Smoothening magnetization distribution\n");
00545
00546
00547 Vec Vtmp;
00548 ierr = VecDuplicate(gdata->M,&Vtmp);CHKERRQ(ierr);
00549
00550 for (int i=0;i<NSMOOTHSTEPS;i++) {
00551
00552 ierr = FieldSmoothApply(gdata,gdata->M,Vtmp);
00553 ierr = FieldSmoothApply(gdata,Vtmp,gdata->M);
00554
00555
00556 PetscReal devNorm;
00557 devNorm=RenormVec(gdata->M,1.0,0.0,ND);
00558 }
00559
00560 ierr = VecDestroy(Vtmp);CHKERRQ(ierr);
00561 }
00562
00563 if (gdata->mode!=53) {
00564 ierr = Cart2Sphere(gdata->M,Msphere);CHKERRQ(ierr);
00565 }
00566
00567 if (NSMOOTHSTEPS>0) {
00568 PetscPrintf(MPI_COMM_WORLD,"Restarting energy minimization with TAO after smoothening...\n");
00569 ierr = TaoSolveApplication(taoapp,tao);CHKERRQ(ierr);
00570 }
00571 #endif
00572
00573 if (gdata->mode!=53) {
00574 Vec sol;
00575 ierr = TaoAppGetSolutionVec(taoapp,&sol);CHKERRQ(ierr);
00576 assert(sol==Msphere);
00577 }
00578
00579 int iter;
00580 double ff,gnorm,cnorm,xdiff;
00581 TaoTerminateReason reason;
00582 ierr = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,&cnorm,&xdiff,&reason);CHKERRQ(ierr);
00583 if (reason <= 0 ){
00584 PetscPrintf(MPI_COMM_WORLD,"TAO solver did NOT converge - ");
00585
00586
00587
00588 gdata->equil=0;
00589 }
00590 else {
00591 PetscPrintf(MPI_COMM_WORLD,"TAO solver converged - ");
00592 gdata->equil++;
00593 }
00594 PetscPrintf(MPI_COMM_WORLD,
00595 "TerminationReason: %i\n"
00596 "it: %i, f: %4.2e, res: %4.2e, cnorm: %4.2e, step %4.2e\n",
00597 reason,
00598 iter,ff,gnorm,cnorm,xdiff
00599 );
00600
00601 ierr = TaoView(tao);CHKERRQ(ierr);
00602
00603
00604 ierr = Htot(gdata);CHKERRQ(ierr);
00605 ierr = Htot_Energy(gdata);CHKERRQ(ierr);
00606
00607 MagparFunctionLogReturn(0);
00608 }
00609