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: mytscreatepvode.c 2681 2009-07-31 04:30:53Z 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;
00038 N_Vector nvu;
00039 realtype abstol;
00040 realtype reltol;
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;
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
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
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
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
00215 SETERRQ(PETSC_ERR_SUP,"Not yet implemented!\n");
00216 break;
00217 default:
00218
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
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
00257 RenormVec(gdata->M,1.0,D_EPS,ND);
00258
00259 #else
00260
00261
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
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
00352
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
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
00375 ierr = PVodeInit(gdata);CHKERRQ(ierr);
00376
00377 MagparFunctionLogReturn(0);
00378 }
00379