eminisolve.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: 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 /* undef DAMP to use conventional gradient
00043    define DAMP to use damping term for gradient/line search
00044 */
00045 #undef DAMP
00046 
00047 /* TODO: problems with using Cartesian coordinates in micromagnetic model:
00048    |M|==1 is not preserved! M "runs away"! - need penalty term or such
00049    works great with nonlinear model
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     /* FIXME: try to avoid copying vectors here
00084        can we do gdata->M=X instead?
00085        same for gdata->VHtot=G ?
00086        search for all VecCopy and try to do something smarter
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   /* use damping term instead of field for gradient/line search */
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     /* skip non-magnetic vertices */
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     /* use small parameter to keep residual for TAO small */
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 /* Test gradient with these options:
00146   -tao_method tao_fd_test
00147   -tao_test_gradient
00148   -tao_test_display
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   /* DEBUG
00180   PetscPrintf(PETSC_COMM_WORLD,"deb01 X\n");
00181   VecView(X,PETSC_VIEWER_STDOUT_WORLD);
00182   PetscPrintf(PETSC_COMM_WORLD,"deb01 VHtot\n");
00183   VecView(gdata->VHtot,PETSC_VIEWER_STDOUT_WORLD);
00184   PetscPrintf(PETSC_COMM_WORLD,"deb01 G\n");
00185   VecView(G,PETSC_VIEWER_STDOUT_WORLD);
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   /* test based on torque dM/dt - does not work well
00208      torque>0 in energy minimum!?
00209      however torque->0 is not (strict) energy minimum
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     /* count how long we have been in equilibrium in succession */
00231     equil++;
00232   }
00233   else {
00234     ierr = TaoSetTerminationReason(tao,TAO_CONTINUE_ITERATING);CHKERRQ(ierr);
00235     /* reset counter if we are not in equilibrium */
00236     equil=0;
00237   }
00238 
00239 #else
00240   /* test based on change in energy */
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   /* Etotold-gdata->Etot < 0 when energy increases
00273      This can happen around energy minimum, due to numerical errors.
00274      We allow it as long as |delta E|<tol.
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       /* count how long we have been in equilibrium in succession */
00285       equil++;
00286     }
00287   }
00288   else {
00289     ierr = TaoSetTerminationReason(tao,TAO_CONTINUE_ITERATING);CHKERRQ(ierr);
00290     /* reset counter if we are not in equilibrium */
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   /* set default for tao_lmm_vectors */
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       /* if the magnetic moment of the material is 0 */
00348       if (ms==0.0) {
00349         /* set magnetization to 0 */
00350         /* This is necessary, because we set Hnonlin=0 for Ms=0 or mu=1 / chi=0 in hnonlin.c */
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   PetscPrintf(PETSC_COMM_WORLD,"Activating TAO Monitor\n");
00378   ierr = PetscOptionsSetValue("-tao_monitor","");CHKERRQ(ierr);
00379 */
00380 
00381   /* set defaults for line search options */
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   /* Check for TAO command line options (listed here for consistency check with allopt.txt) */
00405   /* evaluates
00406       PetscOptions "-tao_fatol"
00407       PetscOptions "-tao_frtol"
00408   */
00409 /* defined in tao-1.9/src/petsctao/include/taopetsc.h - may not be supported in future releases
00410   ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);
00411 */
00412   ierr = TaoSetOptions(taoapp,tao);CHKERRQ(ierr);
00413 
00414   /*
00415     disturb magnetization distribution with small random vector;
00416     larger distortions require many more interations in TaoSolve
00417     and do not help, because the solution is less accurate
00418     (relative tolerance!)
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   /* To View TAO solver information use */
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     /* distort by tilting magnetization in direction of MxMxH */
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       /* TODO: can lead to Inf/NaNs if "damping constant" is too large */
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     /* randomly distort M */
00491     ierr = DistortVec(rctx,gdata->M,magdist);CHKERRQ(ierr);
00492 
00493     /* renormalize M */
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     PetscReal *ta_X, *ta_VMs3;
00510     ierr = VecGetArray(Msphere,&ta_X);CHKERRQ(ierr);
00511     ierr = VecGetArray(gdata->VMs3,&ta_VMs3);CHKERRQ(ierr);
00512     for (int i=0; i<gdata->ln_vert; i++) {
00513       ta_X[2*i+0] *= ta_VMs3[ND*i];
00514       ta_X[2*i+1] *= ta_VMs3[ND*i];
00515     }
00516     ierr = VecRestoreArray(Msphere,&ta_X);CHKERRQ(ierr);
00517     ierr = VecRestoreArray(gdata->VMs3,&ta_VMs3);CHKERRQ(ierr);
00518 */
00519 #endif
00520   }
00521 
00522   ierr = TaoSolveApplication(taoapp,tao);CHKERRQ(ierr);
00523 
00524 #ifdef ADDONS
00525 
00526 #ifdef DEBUG
00527   /* recalculate fields and energies */
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     /* smoothen magnetization distribution */
00547     Vec Vtmp;
00548     ierr = VecDuplicate(gdata->M,&Vtmp);CHKERRQ(ierr);
00549 
00550     for (int i=0;i<NSMOOTHSTEPS;i++) {
00551       /* perform smoothening (double) step */
00552       ierr = FieldSmoothApply(gdata,gdata->M,Vtmp);
00553       ierr = FieldSmoothApply(gdata,Vtmp,gdata->M);
00554 
00555       /* renormalize M (distorted due to smoothening) */
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   /* store this info in log file!? */
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     /* continue even if TAO fails
00586        SETERRQ(PETSC_ERR_LIB,"TAO failed! Exiting!\n");
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   /* otherwise Htot=0 */
00604   ierr = Htot(gdata);CHKERRQ(ierr);
00605   ierr = Htot_Energy(gdata);CHKERRQ(ierr);
00606 
00607   MagparFunctionLogReturn(0);
00608 }
00609 

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