precond.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: 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   /* get ownership range */
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       /* we could also use MPI_Send/Recv to get everything just to the first processor */
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   /* get ownership range */
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   /* calculate number of nonzeroes for each line */
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     /* add 0 to diagonal matrix elements to enforce their allocation
00154        otherwise there might occur problems with MatScale/MatShift in Precond
00155        due to DIFFERENT_NONZERO_PATTERN
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   /* compute the Jacobian */
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   /* set default options */
00185 
00186   /* evaluates
00187       PetscOptions "-psolve_ksp_type"
00188       PetscOptions "-psolve_ksp_atol"
00189       PetscOptions "-psolve_ksp_divtol"
00190       PetscOptions "-psolve_ksp_rtol"
00191                    "-psolve_pc_type"
00192                    "-psolve_sub_pc_type"
00193                    "-psolve_sub_pc_factor_shift_positive_definite"
00194                    "-psolve_pc_type"
00195                    "-psolve_pc_factor_shift_positive_definite"
00196   */
00197 
00198   PetscTruth flg;
00199   char str[256];
00200   /* by default use conjugate gradient solver */
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     /* by default use BlockJacobi preconditioner */
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     /* by default use ILU preconditioner on subblocks */
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     /* by default use Manteuffel shift */
00259 /*
00260 #undef ostr
00261 #define ostr "-psolve_sub_pc_factor_shift_positive_definite"
00262     ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
00263     if (flg!=PETSC_TRUE) {
00264       sprintf(str,"1");
00265       PetscOptionsSetValue(ostr,str);
00266       PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
00267     }
00268 */
00269   }
00270   else {
00271     /* by default use ILU preconditioner */
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     /* by default use Manteuffel shift */
00282 /*
00283 #undef ostr
00284 #define ostr "-psolve_pc_factor_shift_positive_definite"
00285     ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
00286     if (flg!=PETSC_TRUE) {
00287       sprintf(str,"1");
00288       PetscOptionsSetValue(ostr,str);
00289       PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
00290     }
00291 */
00292   }
00293 
00294   /* use prefix for options */
00295   ierr = KSPSetOptionsPrefix(ksp_Jac,"psolve_");CHKERRQ(ierr);
00296 
00297   /* possibly override defaults with settings from user defined options */
00298   ierr = KSPSetFromOptions(ksp_Jac);CHKERRQ(ierr);
00299 
00300 /*
00301   ierr = KSPSetInitialGuessKnoll(ksp_Jac,PETSC_TRUE);CHKERRQ(ierr);
00302   ierr = KSPSetInitialGuessNonzero(ksp_Jac,PETSC_TRUE);CHKERRQ(ierr);
00303 */
00304 
00305   /* KSPSetUp only necessary to get more useful output from KSPView */
00306   /* FIXME: bug(?) in PETSc 2.2.1: for size==1 leads to error "Zero pivot row 0" */
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   /* jok - TRUE means reuse current Jacobian else recompute Jacobian */
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     /* make PETSc vector tmpy point to PVODE vector y */
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     /* TODO: make Htot not calculate demag - we don't use it for Jacobian, could save (a little) time */
00348     ierr = Htot(gdata);CHKERRQ(ierr);
00349 
00350     /* compute the Jacobian */
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     /* copy the Jacobian matrix */
00356     ierr = MatCopy(vode_Jac,vode_pmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
00357 
00358     *jcurPtr = TRUE;
00359   }
00360 
00361   /* construct I-gamma*Jac  */
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 /* adapted from $PETSC_DIR/src/ts/impls/implicit/pvode/petscpvode.c */
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   /* as defined in call to CVSpgmr in pvode_init: 1==left, 2==right */
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 /* no difference in performance visible with mumag std prob. #4
00390    between abstol=0.1*delta*abstol and reltol=0.01
00391    because just 1 interation is usually sufficient
00392    (except if the Jacobian has been recalculated)
00393 
00394   KSP ksp_Jac;
00395   ierr = KSPSetTolerances(ksp_Jac,0.01,0.0,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
00396 
00397   PetscInfo3(0,"set atol to: 0.1*delta*ewt[0]=%g*%g=%g\n",
00398     delta,NV_DATA_P(ewt)[0],0.1*delta*NV_DATA_P(ewt)[0]
00399   );
00400   ierr = KSPSetTolerances(ksp_Jac,0.0,0.1*delta*NV_DATA_P(ewt)[0],PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
00401 */
00402 
00403 /* TODO: useful/necessary before every call to KSPSolve?
00404 */
00405 /*
00406   ierr = KSPSetInitialGuessKnoll(ksp_Jac,PETSC_TRUE);CHKERRQ(ierr);
00407   ierr = KSPSetInitialGuessNonzero(ksp_Jac,PETSC_TRUE);CHKERRQ(ierr);
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 }

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