mytscreatepvode.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: mytscreatepvode.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";
00026 
00027 #include "llg.h"
00028 #include "util/util.h"
00029 
00030 #if SUNDIALS_VERSION >= 230
00031 #include "cvode/cvode_diag.h"
00032 #include "cvode/cvode_spgmr.h"
00033 #else
00034 #error "SUNDIALS versions lower/older than 2.3.0 are not supported"
00035 #endif
00036 
00037 void     *cvode_mem=PETSC_NULL; /* data structures for PVODE */
00038 N_Vector nvu;        /* vector for independent variable (i.e. magnetization) */
00039 realtype abstol;     /* absolute tolerance */
00040 realtype reltol;     /* relative tolerance */
00041 
00042 int PVodeInit(GridData *gdata)
00043 {
00044   int          flag;
00045   PetscTruth   flg;
00046 
00047   MagparFunctionLogBegin;
00048 
00049   char str[256];
00050   ierr = PetscOptionsGetString(PETSC_NULL,"-ts_pvode_type",str,256,&flg);CHKERRQ(ierr);
00051   if (flg!=PETSC_TRUE) {
00052     sprintf(str,"bdf");
00053     PetscPrintf(PETSC_COMM_WORLD,
00054       "Option -ts_pvode_type not found, using default value: %s\n",
00055       str
00056     );
00057   }
00058 
00059   int lmmn=-1;
00060   ierr = PetscStrcmp(str,"adams",&flg);CHKERRQ(ierr);
00061   if (flg) lmmn=CV_ADAMS;
00062   ierr = PetscStrcmp(str,"bdf",&flg);CHKERRQ(ierr);
00063   if (flg) lmmn=CV_BDF;
00064   if (lmmn != CV_ADAMS && lmmn != CV_BDF) {
00065     SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Unknown time integration method.\n");
00066   }
00067   cvode_mem = CVodeCreate(lmmn,lmmn==CV_ADAMS ? CV_FUNCTIONAL : CV_NEWTON);
00068   PetscPrintf(PETSC_COMM_WORLD,"PVODE time integration method: %s\n",str);
00069 
00070   int maxorder;
00071   ierr = PetscOptionsGetInt(PETSC_NULL,"-maxorder",(PetscInt*)&maxorder,&flg);CHKERRQ(ierr);
00072   if (flg!=PETSC_TRUE) {
00073     maxorder=2;
00074     PetscPrintf(PETSC_COMM_WORLD,
00075       "Option -maxorder not found, using default value: %i\n",
00076       maxorder
00077     );
00078   }
00079   flag=CVodeSetMaxOrd(cvode_mem,maxorder);
00080 
00081   PetscReal delt; /* ratio between tolerance of nonlinear and linear iterations */
00082   ierr = PetscOptionsGetReal(PETSC_NULL,"-ts_pvode_linear_tolerance",&delt,&flg);CHKERRQ(ierr);
00083   if (flg!=PETSC_TRUE) {
00084     delt=0.05;
00085     PetscPrintf(PETSC_COMM_WORLD,
00086       "Option -ts_pvode_linear_tolerance not found, using default value: %g\n",
00087       delt
00088     );
00089   }
00090   PetscPrintf(PETSC_COMM_WORLD,"delt: %g\n",delt);
00091 
00092   PetscReal dtime;
00093   ierr = PetscOptionsGetReal(PETSC_NULL,"-ts_dt",&dtime,&flg);CHKERRQ(ierr);
00094   if (flg!=PETSC_TRUE) {
00095     dtime = gdata->tscale*1e9/100;
00096     PetscPrintf(PETSC_COMM_WORLD,
00097       "Option -ts_dt not found, using default value: %g\n",
00098       dtime
00099     );
00100   }
00101 
00102   PetscReal dtmin;
00103   ierr = PetscOptionsGetReal(PETSC_NULL,"-mintimestep",&dtmin,&flg);CHKERRQ(ierr);
00104   if (flg!=PETSC_TRUE) {
00105     dtmin=0.0;
00106     PetscPrintf(PETSC_COMM_WORLD,
00107       "Option -mintimestep not found, using default value: %g\n",
00108       dtmin
00109     );
00110   }
00111 
00112   PetscReal dtmax;
00113   ierr = PetscOptionsGetReal(PETSC_NULL,"-maxtimestep",&dtmax,&flg);CHKERRQ(ierr);
00114   if (flg!=PETSC_TRUE) {
00115     dtmax=1;
00116     PetscPrintf(PETSC_COMM_WORLD,
00117       "Option -maxtimestep not found, using default value: %g\n",
00118       dtmax
00119     );
00120   }
00121 
00122   dtime /= 1e9*gdata->tscale;
00123   dtmin /= 1e9*gdata->tscale;
00124   dtmax /= 1e9*gdata->tscale;
00125 
00126   /* dtime must not be zero! otherwise PVODE does not start properly */
00127   dtime = dtime < D_EPS ? dtmin : dtime;
00128   dtmin = dtmin < 0 ? dtime : dtmin;
00129   dtmax = dtmax < 1000 ? dtmax : 1000;
00130 
00131   PetscPrintf(PETSC_COMM_WORLD,"initial time step: %g ns\n",dtime*gdata->tscale*1e9);
00132   PetscPrintf(PETSC_COMM_WORLD,"minimum time step: %g ns\n",dtmin*gdata->tscale*1e9);
00133   PetscPrintf(PETSC_COMM_WORLD,"maximum time step: %g ns\n",dtmax*gdata->tscale*1e9);
00134   PetscPrintf(PETSC_COMM_WORLD,"reduced initial time step: %g (1)\n",dtime);
00135   PetscPrintf(PETSC_COMM_WORLD,"reduced minimum time step: %g (1)\n",dtmin);
00136   PetscPrintf(PETSC_COMM_WORLD,"reduced maximum time step: %g (1)\n",dtmax);
00137 
00138   flag=CVodeSetInitStep(cvode_mem,dtime);
00139   flag=CVodeSetMinStep(cvode_mem,dtmin);
00140   flag=CVodeSetMaxStep(cvode_mem,dtmax);
00141 
00142   CVodeSetFdata(cvode_mem,gdata);
00143 
00144   flag = CVodeMalloc(
00145     cvode_mem,
00146     RHSfunction,
00147     gdata->time,
00148     nvu,
00149     CV_SS,
00150     reltol,
00151     &abstol
00152   );
00153   if (flag != CV_SUCCESS) {
00154     SETERRQ1(PETSC_ERR_LIB,"CVodeMalloc failed: %i\n",flag);
00155   }
00156 
00157   int precon;
00158   ierr = PetscOptionsGetInt(PETSC_NULL,"-precon",(PetscInt*)&precon,&flg);CHKERRQ(ierr);
00159   if (flg!=PETSC_TRUE) {
00160     precon=1;
00161     PetscPrintf(PETSC_COMM_WORLD,
00162       "Option -precon not found, using default value: %i\n",
00163       precon
00164     );
00165   }
00166   PetscPrintf(PETSC_COMM_WORLD,"preconditioning method: %i\n",precon);
00167 
00168   switch (precon) {
00169     case 0:
00170       PetscPrintf(PETSC_COMM_WORLD,"preconditioning off\n");
00171       flag = CVDiag(cvode_mem);
00172     break;
00173     case 1:
00174       /* Jacobi preconditioning - use Jtimes only, if demag is zero */
00175       PetscPrintf(PETSC_COMM_WORLD,"Jacobi preconditioning on\n");
00176 
00177       int maxl;        
00178       ierr = PetscOptionsGetInt(PETSC_NULL,"-maxl",(PetscInt*)&maxl,&flg);CHKERRQ(ierr);
00179       if (flg!=PETSC_TRUE) {
00180         maxl=300;
00181         PetscPrintf(PETSC_COMM_WORLD,
00182           "Option -maxl not found, using default value: %i\n",
00183           maxl
00184         );
00185       }
00186       PetscPrintf(PETSC_COMM_WORLD,"maxl: %i\n",maxl);
00187       flag = CVSpgmr(cvode_mem,PREC_LEFT,(int)maxl);
00188 
00189       /* initialize Jacobi preconditioning */
00190       ierr = Precond_Init(gdata);CHKERRQ(ierr);
00191 
00192 #if SUNDIALS_VERSION >= 230
00193       flag = CVSpilsSetPrecType(cvode_mem,PREC_LEFT);
00194       flag = CVSpilsSetGSType(cvode_mem,MODIFIED_GS);
00195       flag = CVSpilsSetPreconditioner(cvode_mem,Precond,PSolve,gdata);
00196 #else
00197 #error "SUNDIALS versions lower/older than 2.3.0 are not supported"
00198 #endif
00199       if (gdata->VHdem==PETSC_NULL) {
00200 #if SUNDIALS_VERSION >= 230
00201         flag = CVSpilsSetJacTimesVecFn(cvode_mem,Jtimes,gdata);
00202 #else
00203 #error "SUNDIALS versions lower/older than 2.3.0 are not supported"
00204 #endif
00205         PetscPrintf(PETSC_COMM_WORLD,"Hdemag off: PVODE uses Jtimes function\n");
00206       }
00207       else {
00208         PetscPrintf(PETSC_COMM_WORLD,
00209           "Hdemag on: PVODE uses internal difference quotient approximation for Jacobian*vector\n"
00210         );
00211       }
00212     break;
00213     case 2:
00214       /* BBD preconditioning */
00215       SETERRQ(PETSC_ERR_SUP,"Not yet implemented!\n");
00216     break;
00217     default:
00218       /* no preconditioning */
00219       PetscPrintf(PETSC_COMM_WORLD,
00220         "PVODE uses internal difference quotient approximation for Jacobian*vector\n"
00221       );
00222       flag = CVDiag(cvode_mem);
00223     break;
00224   }
00225 
00226   if (flag == CV_SUCCESS) {
00227     PetscPrintf(PETSC_COMM_WORLD,"PVODE successfully initialized\n");
00228   }
00229   else {
00230     SETERRQ(PETSC_ERR_MEM,"CVSpgmr failed.\n");
00231   }
00232 
00233   MagparFunctionLogReturn(0);
00234 }
00235 
00236 
00237 int PVodeReInit(GridData *gdata)
00238 {
00239   MagparFunctionInfoBegin;
00240 
00241 #if 1
00242 
00243   /* TODO: remove this section between VecGet/RestoreArray after sufficient testing */
00244 
00245   PetscReal *tmp;
00246   ierr = VecGetArray(gdata->M,&tmp);CHKERRQ(ierr);
00247 
00248 #ifdef UNIPROC
00249   assert(tmp==NV_DATA_S(nvu));
00250 #else
00251   assert(tmp==NV_DATA_P(nvu));
00252 #endif
00253 
00254   ierr = VecRestoreArray(gdata->M,&tmp);CHKERRQ(ierr);
00255 
00256   /* renormalize M since we just reinitialized PVODE */
00257   RenormVec(gdata->M,1.0,D_EPS,ND);
00258 
00259 #else
00260 
00261 /* TODO: remove this section after sufficient testing */
00262 
00263 #ifdef UNIPROC
00264   VecPlaceArray(gdata->M,NV_DATA_S(nvu));
00265 #else
00266   VecPlaceArray(gdata->M,NV_DATA_P(nvu));
00267 #endif
00268 
00269   /* renormalize M since we just reinitialized PVODE */
00270   RenormVec(gdata->M,1.0,D_EPS,ND);
00271 
00272   ierr=VecResetArray(gdata->M);
00273 #endif
00274 
00275   PetscInfo1(0,"t=%g ns: renormalized M\n",gdata->time*gdata->tscale*1e9);
00276 
00277   int flag;
00278 
00279 #undef RESETDTIME
00280 #ifdef RESETDTIME
00281   realtype dtime;
00282   flag=CVodeGetCurrentStep(cvode_mem,&dtime);
00283 #endif
00284 
00285   flag =
00286     CVodeReInit(
00287       cvode_mem,
00288       RHSfunction,
00289       gdata->time,
00290       nvu,
00291       CV_SS,
00292       reltol,
00293       &abstol
00294     );
00295   if (flag != CV_SUCCESS) {
00296     SETERRQ1(PETSC_ERR_LIB,"CVodeReInit failed: %i\n", flag);
00297   }
00298   PetscInfo(0,"PVODE successfully reinitialized\n");
00299 
00300 #ifdef RESETDTIME
00301   flag=CVodeSetInitStep(cvode_mem,dtime);
00302 #endif
00303 
00304   MagparFunctionProfReturn(0);
00305 }
00306 
00307 
00308 int myTSCreatePVode(GridData *gdata)
00309 {
00310   PetscTruth flg;
00311 
00312   MagparFunctionLogBegin;
00313 
00314   int rank,size;
00315   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00316   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00317 
00318   ierr = PetscOptionsGetReal(PETSC_NULL,"-ts_pvode_rtol",&reltol,  &flg);CHKERRQ(ierr);
00319   if (flg!=PETSC_TRUE) {
00320     reltol=1e-6;
00321     PetscPrintf(PETSC_COMM_WORLD,
00322       "Option -ts_pvode_rtol not found, using default value: %g\n",
00323       reltol
00324     );
00325   }
00326   PetscPrintf(PETSC_COMM_WORLD,"reltol:      %g\n",reltol);
00327 
00328   ierr = PetscOptionsGetReal(PETSC_NULL,"-ts_pvode_atol",&abstol,  &flg);CHKERRQ(ierr);
00329   if (flg!=PETSC_TRUE) {
00330     abstol=1e-6;
00331     PetscPrintf(PETSC_COMM_WORLD,
00332       "Option -ts_pvode_atol not found, using default value: %g\n",
00333       abstol
00334     );
00335   }
00336   PetscPrintf(PETSC_COMM_WORLD,"abstol:      %g\n",abstol);
00337 
00338   int neq, local_N;
00339   ierr = VecGetSize(gdata->M,(PetscInt*)&neq);CHKERRQ(ierr);
00340   ierr = VecGetLocalSize(gdata->M,(PetscInt*)&local_N);CHKERRQ(ierr);
00341 
00342   assert(neq==ND*gdata->n_vert);
00343   assert(local_N==ND*gdata->ln_vert);
00344 
00345 #ifdef UNIPROC
00346   nvu=N_VNew_Serial(neq);
00347 #else
00348   nvu=N_VNew_Parallel(PETSC_COMM_WORLD,local_N,neq);
00349 #endif
00350 
00351   /* create new vector with PVODE's array */
00352   /* real data must be of same type as PETSc data (usually double) !!! */
00353   Vec  Mnew;
00354   ierr = VecCreateMPIWithArray(
00355     PETSC_COMM_WORLD,
00356     local_N,
00357     neq,
00358 #ifdef UNIPROC
00359     NV_DATA_S(nvu),
00360 #else
00361     NV_DATA_P(nvu),
00362 #endif
00363     &Mnew
00364   );CHKERRQ(ierr);
00365 
00366   /* copy current magnetization distribution into the new vector */
00367   VecCopy(gdata->M,Mnew);
00368   VecDestroy(gdata->M);
00369   gdata->M=Mnew;
00370   ierr = VecSetBlockSize(gdata->M,ND);CHKERRQ(ierr);
00371 
00372   cvode_mem=NULL;
00373 
00374   /* initialize PVode */
00375   ierr = PVodeInit(gdata);CHKERRQ(ierr);
00376 
00377   MagparFunctionLogReturn(0);
00378 }
00379 

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