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: precond.c 3013 2010-03-26 16:12:31Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027 #include "llg.h"
00028 #include "field/field.h"
00029 #include "util/util.h"
00030
00031 static KSP ksp_Jac;
00032 static Mat vode_Jac;
00033 static Vec vode_w1,vode_w2;
00034 static Mat vode_pmat;
00036 int Precond_Init(GridData *gdata)
00037 {
00038 MagparFunctionLogBegin;
00039
00040 int rank,size;
00041 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00042 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00043
00044 PetscPrintf(PETSC_COMM_WORLD,"Calculating sparsity pattern for matrix preallocation\n");
00045
00046
00047 int low,high;
00048 ierr = VecGetOwnershipRange(gdata->elevol,&low,&high);CHKERRQ(ierr);
00049
00050 int *elevertall;
00051
00052 #ifdef HAVE_ELEVERTALL
00053 elevertall=gdata->elevertall;
00054 #else
00055 if (size>1) {
00056 ierr = PetscMalloc(NV*gdata->n_ele*sizeof(int),&elevertall);CHKERRQ(ierr);
00057 ierr = PetscMemcpy(elevertall+NV*low,gdata->elevert,NV*gdata->ln_ele*sizeof(int));CHKERRQ(ierr);
00058
00059 int *p;
00060 p=elevertall;
00061
00062 int recvcount[size];
00063 recvcount[rank]=NV*gdata->ln_ele;
00064 for (int i=0;i<size;i++) {
00065 ierr = MPI_Bcast(recvcount+i,1,MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
00066
00067 ierr = MPI_Bcast(p,recvcount[i],MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
00068 p+=recvcount[i];
00069 }
00070 assert(p==elevertall+NV*gdata->n_ele);
00071 }
00072 else {
00073 elevertall=gdata->elevert;
00074 }
00075 #endif
00076
00077 int *ia,*ja;
00078 #ifdef HAVE_M2IJ
00079 ia=gdata->mesh2nodal_ia;
00080 ja=gdata->mesh2nodal_ja;
00081 #else
00082 ierr = Mesh2Nodal(gdata->n_ele,gdata->n_vert,elevertall,&ia,&ja);CHKERRQ(ierr);
00083 #endif
00084
00085 #ifndef HAVE_ELEVERTALL
00086 if (size>1) {
00087 ierr = PetscFree(elevertall);CHKERRQ(ierr);
00088 }
00089 #endif
00090
00091
00092 ierr = VecGetOwnershipRange(gdata->M,&low,&high);CHKERRQ(ierr);
00093 low/=ND;
00094 high/=ND;
00095
00096 int *d_nnz;
00097 int *o_nnz;
00098 ierr = PetscMalloc(ND*(high-low)*sizeof(int),&d_nnz);CHKERRQ(ierr);
00099 ierr = PetscMalloc(ND*(high-low)*sizeof(int),&o_nnz);CHKERRQ(ierr);
00100
00101 for (int i=low;i<high;i++) {
00102 for (int j=0;j<ND;j++) {
00103 d_nnz[ND*(i-low)+j]=ND;
00104 o_nnz[ND*(i-low)+j]=0;
00105 for (int l=ia[i];l<ia[i+1];l++) {
00106 if (ja[l]>=low && ja[l]<high) {
00107 d_nnz[ND*(i-low)+j]+=ND;
00108 }
00109 else {
00110 o_nnz[ND*(i-low)+j]+=ND;
00111 }
00112 }
00113 }
00114 }
00115
00116 #if 0
00117 #define PREALLOC_DG 25
00118 #define PREALLOC_OD 15
00119 PetscPrintf(PETSC_COMM_WORLD,"Creating matrix with fixed preallocation\n");
00120
00121 ierr = MatCreateMPIAIJ(
00122 PETSC_COMM_WORLD,
00123 ND*gdata->ln_vert,ND*gdata->ln_vert,
00124 ND*gdata->n_vert,ND*gdata->n_vert,
00125 PREALLOC_DG,0,PREALLOC_OD,0,
00126 &vode_Jac
00127 );CHKERRQ(ierr);
00128 #else
00129 PetscPrintf(PETSC_COMM_WORLD,"Creating matrix with calculated preallocation\n");
00130 ierr = MatCreateMPIAIJ(
00131 PETSC_COMM_WORLD,
00132 ND*gdata->ln_vert,ND*gdata->ln_vert,
00133 ND*gdata->n_vert,ND*gdata->n_vert,
00134 0,d_nnz,0,o_nnz,
00135 &vode_Jac
00136 );CHKERRQ(ierr);
00137 #endif
00138 ierr = MatSetFromOptions(vode_Jac);CHKERRQ(ierr);
00139
00140 #ifndef HAVE_M2IJ
00141 ierr = PetscFree(ia);CHKERRQ(ierr);
00142 ierr = PetscFree(ja);CHKERRQ(ierr);
00143 #endif
00144 ierr = PetscFree(d_nnz);CHKERRQ(ierr);
00145 ierr = PetscFree(o_nnz);CHKERRQ(ierr);
00146
00147 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,ND*gdata->ln_vert,ND*gdata->n_vert,PETSC_NULL,&vode_w1);CHKERRQ(ierr);
00148 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,ND*gdata->ln_vert,ND*gdata->n_vert,PETSC_NULL,&vode_w2);CHKERRQ(ierr);
00149
00150 PetscPrintf(PETSC_COMM_WORLD,"Calculating matrix elements...\n");
00151
00152 for (int i=0;i<gdata->ln_vert;i++) {
00153
00154
00155
00156
00157
00158 int rows[ND];
00159 rows[0]=gdata->vertl2g[i]*ND;
00160 rows[1]=rows[0]+1;
00161 rows[2]=rows[0]+2;
00162
00163 ierr = MatSetValue(vode_Jac,rows[0],rows[0],0.0,ADD_VALUES);CHKERRQ(ierr);
00164 ierr = MatSetValue(vode_Jac,rows[1],rows[1],0.0,ADD_VALUES);CHKERRQ(ierr);
00165 ierr = MatSetValue(vode_Jac,rows[2],rows[2],0.0,ADD_VALUES);CHKERRQ(ierr);
00166 }
00167
00168
00169 MatStructure str2 = DIFFERENT_NONZERO_PATTERN;
00170 ierr = myLLGJacobian(PETSC_NULL,gdata->time,gdata->M,&vode_Jac,&vode_Jac,&str2,gdata);CHKERRQ(ierr);
00171
00172 ierr = MatDuplicate(vode_Jac,MAT_DO_NOT_COPY_VALUES,&vode_pmat);CHKERRQ(ierr);
00173 #ifdef SUPERLUX
00174 ierr = MatConvert(vode_Jac,MATSUPERLU_DIST,&vode_Jac);CHKERRQ(ierr);
00175 ierr = MatConvert(vode_pmat,MATSUPERLU_DIST,&vode_pmat);CHKERRQ(ierr);
00176 ierr = MatDestroy(vode_Jac);CHKERRQ(ierr);
00177 #endif
00178
00179 PetscPrintf(PETSC_COMM_WORLD,"Initializing linear solver...\n");
00180
00181 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp_Jac);CHKERRQ(ierr);
00182 ierr = KSPSetOperators(ksp_Jac,vode_Jac,vode_Jac,SAME_PRECONDITIONER);CHKERRQ(ierr);
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198 PetscTruth flg;
00199 char str[256];
00200
00201 #undef ostr
00202 #define ostr "-psolve_ksp_type"
00203 ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
00204 if (flg!=PETSC_TRUE) {
00205 sprintf(str,"gmres");
00206 PetscOptionsSetValue(ostr,str);
00207 PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
00208 }
00209
00210 #undef ostr
00211 #define ostr "-psolve_ksp_atol"
00212 ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
00213 if (flg!=PETSC_TRUE) {
00214 sprintf(str,"1e-7");
00215 PetscOptionsSetValue(ostr,str);
00216 PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
00217 }
00218
00219 #undef ostr
00220 #define ostr "-psolve_ksp_rtol"
00221 ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
00222 if (flg!=PETSC_TRUE) {
00223 sprintf(str,"0.01");
00224 PetscOptionsSetValue(ostr,str);
00225 PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
00226 }
00227
00228 #undef ostr
00229 #define ostr "-psolve_ksp_divtol"
00230 ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
00231 if (flg!=PETSC_TRUE) {
00232 sprintf(str,"100");
00233 PetscOptionsSetValue(ostr,str);
00234 PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
00235 }
00236
00237 if (size>1) {
00238
00239 #undef ostr
00240 #define ostr "-psolve_pc_type"
00241 ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
00242 if (flg!=PETSC_TRUE) {
00243 sprintf(str,"bjacobi");
00244 PetscOptionsSetValue(ostr,str);
00245 PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
00246 }
00247
00248
00249 #undef ostr
00250 #define ostr "-psolve_sub_pc_type"
00251 ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
00252 if (flg!=PETSC_TRUE) {
00253 sprintf(str,"ilu");
00254 PetscOptionsSetValue(ostr,str);
00255 PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
00256 }
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269 }
00270 else {
00271
00272 #undef ostr
00273 #define ostr "-psolve_pc_type"
00274 ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
00275 if (flg!=PETSC_TRUE) {
00276 sprintf(str,"ilu");
00277 PetscOptionsSetValue(ostr,str);
00278 PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
00279 }
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292 }
00293
00294
00295 ierr = KSPSetOptionsPrefix(ksp_Jac,"psolve_");CHKERRQ(ierr);
00296
00297
00298 ierr = KSPSetFromOptions(ksp_Jac);CHKERRQ(ierr);
00299
00300
00301
00302
00303
00304
00305
00306
00307 if (size>1) {
00308 ierr = KSPSetUp(ksp_Jac);CHKERRQ(ierr);
00309 }
00310
00311 PetscPrintf(PETSC_COMM_WORLD,"KSP for preconditioning: P*z=r:\n");
00312 #if defined(SUPERLU)
00313 PetscReal rtol,atol,divtol;
00314 int maxits;
00315 ierr = KSPGetTolerances(ksp_Jac,&rtol,&atol,&divtol,&maxits);CHKERRQ(ierr);
00316 PetscPrintf(PETSC_COMM_WORLD,"rtol: %g atol: %g divtol: %g maxits: %i\n",
00317 rtol,atol,divtol,maxits
00318 );
00319 #else
00320 ierr = KSPView(ksp_Jac,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
00321 #endif
00322
00323 MagparFunctionLogReturn(0);
00324 }
00325
00326
00327 int Precond(realtype t,N_Vector y,N_Vector fy,booleantype jok,booleantype *jcurPtr,realtype gamma,void *P_data,N_Vector vtemp1,N_Vector vtemp2,N_Vector vtemp3)
00328 {
00329 GridData *gdata = (GridData*)P_data;
00330 MatStructure str = DIFFERENT_NONZERO_PATTERN;
00331
00332 MagparFunctionInfoBegin;
00333
00334
00335 if (jok) {
00336 ierr = MatCopy(vode_pmat,vode_Jac,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
00337 str = SAME_NONZERO_PATTERN;
00338 *jcurPtr = FALSE;
00339 } else {
00340
00341 #ifdef UNIPROC
00342 ierr = VecPlaceArray(gdata->M,NV_DATA_S(y));CHKERRQ(ierr);
00343 #else
00344 ierr = VecPlaceArray(gdata->M,NV_DATA_P(y));CHKERRQ(ierr);
00345 #endif
00346
00347
00348 ierr = Htot(gdata);CHKERRQ(ierr);
00349
00350
00351 ierr = myLLGJacobian(PETSC_NULL,gdata->time,gdata->M,&vode_Jac,&vode_Jac,&str,gdata);CHKERRQ(ierr);
00352
00353 ierr = VecResetArray(gdata->M);CHKERRQ(ierr);
00354
00355
00356 ierr = MatCopy(vode_Jac,vode_pmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
00357
00358 *jcurPtr = TRUE;
00359 }
00360
00361
00362
00363 ierr = MatScale(vode_Jac,-gamma);CHKERRQ(ierr);
00364 ierr = MatShift(vode_Jac,1.0);CHKERRQ(ierr);
00365
00366 ierr = KSPSetOperators(ksp_Jac,vode_Jac,vode_Jac,str);CHKERRQ(ierr);
00367
00368 MagparFunctionProfReturn(0);
00369 }
00370
00371
00372
00373
00374 int PSolve(realtype t,N_Vector y,N_Vector fy,N_Vector r,N_Vector z,realtype gamma,realtype delta,int lr,void *P_data,N_Vector vtemp)
00375 {
00376 MagparFunctionInfoBegin;
00377
00378
00379 assert(lr==1);
00380
00381 #ifdef UNIPROC
00382 ierr = VecPlaceArray(vode_w1,NV_DATA_S(r));CHKERRQ(ierr);
00383 ierr = VecPlaceArray(vode_w2,NV_DATA_S(z));CHKERRQ(ierr);
00384 #else
00385 ierr = VecPlaceArray(vode_w1,NV_DATA_P(r));CHKERRQ(ierr);
00386 ierr = VecPlaceArray(vode_w2,NV_DATA_P(z));CHKERRQ(ierr);
00387 #endif
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409 int its1;
00410 ierr = KSPSolve(ksp_Jac,vode_w1,vode_w2);CHKERRQ(ierr);
00411
00412 ierr = VecResetArray(vode_w1);CHKERRQ(ierr);
00413 ierr = VecResetArray(vode_w2);CHKERRQ(ierr);
00414
00415 ierr = KSPGetIterationNumber(ksp_Jac,(PetscInt*)&its1);CHKERRQ(ierr);
00416 KSPConvergedReason reason;
00417 ierr = KSPGetConvergedReason(ksp_Jac,&reason);CHKERRQ(ierr);
00418
00419 ierr = PetscGetTime(&t_t2);CHKERRQ(ierr);
00420 PetscInfo2(0,"timing: %s took %g s\n",__FUNCT__,t_t2-t_t1);
00421 PetscInfo2(0,"psolve: %i its, reason: %i\n",its1,reason);
00422 if (reason<0) {
00423 PetscReal rnorm;
00424
00425 ierr = KSPGetResidualNorm(ksp_Jac,&rnorm);CHKERRQ(ierr);
00426 PetscPrintf(PETSC_COMM_WORLD,
00427 "Warning: psolve: KSP did not converge - reason: %i, rnorm: %g, its: %i\n",
00428 reason,rnorm,its1
00429 );
00430 }
00431
00432 MagparFunctionProfReturn(0);
00433 }
00434
00435
00436 int Jtimes(N_Vector v,N_Vector Jv,realtype t,N_Vector y,N_Vector fy,void *jac_data,N_Vector tmp)
00437 {
00438 MagparFunctionInfoBegin;
00439
00440 #ifdef UNIPROC
00441 ierr = VecPlaceArray(vode_w1,NV_DATA_S(v));CHKERRQ(ierr);
00442 ierr = VecPlaceArray(vode_w2,NV_DATA_S(Jv));CHKERRQ(ierr);
00443 #else
00444 ierr = VecPlaceArray(vode_w1,NV_DATA_P(v));CHKERRQ(ierr);
00445 ierr = VecPlaceArray(vode_w2,NV_DATA_P(Jv));CHKERRQ(ierr);
00446 #endif
00447
00448 ierr = MatMult(vode_pmat,vode_w1,vode_w2);CHKERRQ(ierr);
00449
00450 ierr = VecResetArray(vode_w1);CHKERRQ(ierr);
00451 ierr = VecResetArray(vode_w2);CHKERRQ(ierr);
00452
00453 MagparFunctionProfReturn(0);
00454 }