precond.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2009 Werner Scholz

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    magpar is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

static char const Id[] = "$Id: precond.c 2699 2009-08-07 18:16:34Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "llg.h"
#include "field/field.h"
#include "util/util.h"

static KSP  ksp_Jac; 
static Mat  vode_Jac; 
static Vec  vode_w1,vode_w2; 
static Mat  vode_pmat; 
int Precond_Init(GridData *gdata)
{
  MagparFunctionLogBegin;

  int rank,size;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  PetscPrintf(PETSC_COMM_WORLD,"Calculating sparsity pattern for matrix preallocation\n");

  /* get ownership range */
  int low,high;
  ierr = VecGetOwnershipRange(gdata->elevol,&low,&high);CHKERRQ(ierr);

  int *elevertall;

#ifdef HAVE_ELEVERTALL
  elevertall=gdata->elevertall;
#else
  if (size>1) {
    ierr = PetscMalloc(NV*gdata->n_ele*sizeof(int),&elevertall);CHKERRQ(ierr);
    ierr = PetscMemcpy(elevertall+NV*low,gdata->elevert,NV*gdata->ln_ele*sizeof(int));CHKERRQ(ierr);

    int *p;
    p=elevertall;

    int recvcount[size];
    recvcount[rank]=NV*gdata->ln_ele;
    for (int i=0;i<size;i++) {
      ierr = MPI_Bcast(recvcount+i,1,MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
      /* we could also use MPI_Send/Recv to get everything just to the first processor */
      ierr = MPI_Bcast(p,recvcount[i],MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
      p+=recvcount[i];
    }
    assert(p==elevertall+NV*gdata->n_ele);
  }
  else {
    elevertall=gdata->elevert;
  }
#endif

  int *ia,*ja;
#ifdef HAVE_M2IJ
  ia=gdata->mesh2nodal_ia;
  ja=gdata->mesh2nodal_ja;
#else
  ierr = Mesh2Nodal(gdata->n_ele,gdata->n_vert,elevertall,&ia,&ja);CHKERRQ(ierr);
#endif

#ifndef HAVE_ELEVERTALL
  if (size>1) {
    ierr = PetscFree(elevertall);CHKERRQ(ierr);
  }
#endif

  /* get ownership range */
  ierr = VecGetOwnershipRange(gdata->M,&low,&high);CHKERRQ(ierr);
  low/=ND;
  high/=ND;

  int *d_nnz;
  int *o_nnz;
  ierr = PetscMalloc(ND*(high-low)*sizeof(int),&d_nnz);CHKERRQ(ierr);
  ierr = PetscMalloc(ND*(high-low)*sizeof(int),&o_nnz);CHKERRQ(ierr);
  /* calculate number of nonzeroes for each line */
  for (int i=low;i<high;i++) {
    for (int j=0;j<ND;j++) {
      d_nnz[ND*(i-low)+j]=ND;
      o_nnz[ND*(i-low)+j]=0;
      for (int l=ia[i];l<ia[i+1];l++) {
        if (ja[l]>=low && ja[l]<high) {
          d_nnz[ND*(i-low)+j]+=ND;
        }
        else {
          o_nnz[ND*(i-low)+j]+=ND;
        }
      }
    }
  }

#if 0
#define PREALLOC_DG 25
#define PREALLOC_OD 15
  PetscPrintf(PETSC_COMM_WORLD,"Creating matrix with fixed preallocation\n");

  ierr = MatCreateMPIAIJ(
    PETSC_COMM_WORLD,
    ND*gdata->ln_vert,ND*gdata->ln_vert,
    ND*gdata->n_vert,ND*gdata->n_vert,
    PREALLOC_DG,0,PREALLOC_OD,0,
    &vode_Jac
  );CHKERRQ(ierr);
#else
  PetscPrintf(PETSC_COMM_WORLD,"Creating matrix with calculated preallocation\n");
  ierr = MatCreateMPIAIJ(
    PETSC_COMM_WORLD,
    ND*gdata->ln_vert,ND*gdata->ln_vert,
    ND*gdata->n_vert,ND*gdata->n_vert,
    0,d_nnz,0,o_nnz,
    &vode_Jac
  );CHKERRQ(ierr);
#endif
  ierr = MatSetFromOptions(vode_Jac);CHKERRQ(ierr);

#ifndef HAVE_M2IJ
  ierr = PetscFree(ia);CHKERRQ(ierr);
  ierr = PetscFree(ja);CHKERRQ(ierr);
#endif
  ierr = PetscFree(d_nnz);CHKERRQ(ierr);
  ierr = PetscFree(o_nnz);CHKERRQ(ierr);

  ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,ND*gdata->ln_vert,ND*gdata->n_vert,PETSC_NULL,&vode_w1);CHKERRQ(ierr);
  ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,ND*gdata->ln_vert,ND*gdata->n_vert,PETSC_NULL,&vode_w2);CHKERRQ(ierr);

  PetscPrintf(PETSC_COMM_WORLD,"Calculating matrix elements...\n");

  for (int i=0;i<gdata->ln_vert;i++) {
    /* add 0 to diagonal matrix elements to enforce their allocation
       otherwise there might occur problems with MatScale/MatShift in Precond
       due to DIFFERENT_NONZERO_PATTERN
    */

    int rows[ND];
    rows[0]=gdata->vertl2g[i]*ND;
    rows[1]=rows[0]+1;
    rows[2]=rows[0]+2;

    ierr = MatSetValue(vode_Jac,rows[0],rows[0],0.0,ADD_VALUES);CHKERRQ(ierr);
    ierr = MatSetValue(vode_Jac,rows[1],rows[1],0.0,ADD_VALUES);CHKERRQ(ierr);
    ierr = MatSetValue(vode_Jac,rows[2],rows[2],0.0,ADD_VALUES);CHKERRQ(ierr);
  }

  /* compute the Jacobian */
  MatStructure str2 = DIFFERENT_NONZERO_PATTERN;
  ierr = myLLGJacobian(PETSC_NULL,gdata->time,gdata->M,&vode_Jac,&vode_Jac,&str2,gdata);CHKERRQ(ierr);

  ierr = MatDuplicate(vode_Jac,MAT_DO_NOT_COPY_VALUES,&vode_pmat);CHKERRQ(ierr);
#ifdef SUPERLUX
  ierr = MatConvert(vode_Jac,MATSUPERLU_DIST,&vode_Jac);CHKERRQ(ierr);
  ierr = MatConvert(vode_pmat,MATSUPERLU_DIST,&vode_pmat);CHKERRQ(ierr);
  ierr = MatDestroy(vode_Jac);CHKERRQ(ierr);
#elif defined(BLOCKSOLVEX) /* takes forever to set up PC */
  ierr = MatConvert(vode_Jac,MATMPIROWBS,&vode_Jac);CHKERRQ(ierr);
  ierr = MatConvert(vode_pmat,MATMPIROWBS,&vode_pmat);CHKERRQ(ierr);
  ierr = MatDestroy(vode_Jac);CHKERRQ(ierr);
#endif

  PetscPrintf(PETSC_COMM_WORLD,"Initializing linear solver...\n");

  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp_Jac);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp_Jac,vode_Jac,vode_Jac,SAME_PRECONDITIONER);CHKERRQ(ierr);

  /* set default options */

  /* evaluates
      PetscOptions "-psolve_ksp_type"
      PetscOptions "-psolve_ksp_atol"
      PetscOptions "-psolve_ksp_divtol"
      PetscOptions "-psolve_ksp_rtol"
                   "-psolve_pc_type"
                   "-psolve_sub_pc_type"
                   "-psolve_sub_pc_factor_shift_positive_definite"
                   "-psolve_pc_type"
                   "-psolve_pc_factor_shift_positive_definite"
  */

  PetscTruth flg;
  char str[256];
  /* by default use conjugate gradient solver */
#undef ostr
#define ostr "-psolve_ksp_type"
  ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    sprintf(str,"gmres");
    PetscOptionsSetValue(ostr,str);
    PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
  }

#undef ostr
#define ostr "-psolve_ksp_atol"
  ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    sprintf(str,"1e-7");
    PetscOptionsSetValue(ostr,str);
    PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
  }

#undef ostr
#define ostr "-psolve_ksp_rtol"
  ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    sprintf(str,"0.01");
    PetscOptionsSetValue(ostr,str);
    PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
  }

#undef ostr
#define ostr "-psolve_ksp_divtol"
  ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    sprintf(str,"100");
    PetscOptionsSetValue(ostr,str);
    PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
  }

  if (size>1) {
    /* by default use BlockJacobi preconditioner */
#undef ostr
#define ostr "-psolve_pc_type"
    ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
    if (flg!=PETSC_TRUE) {
      sprintf(str,"bjacobi");
      PetscOptionsSetValue(ostr,str);
      PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
    }

    /* by default use ILU preconditioner on subblocks */
#undef ostr
#define ostr "-psolve_sub_pc_type"
    ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
    if (flg!=PETSC_TRUE) {
      sprintf(str,"ilu");
      PetscOptionsSetValue(ostr,str);
      PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
    }

    /* by default use Manteuffel shift */
/*
#undef ostr
#define ostr "-psolve_sub_pc_factor_shift_positive_definite"
    ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
    if (flg!=PETSC_TRUE) {
      sprintf(str,"1");
      PetscOptionsSetValue(ostr,str);
      PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
    }
*/
  }
  else {
    /* by default use ILU preconditioner */
#undef ostr
#define ostr "-psolve_pc_type"
    ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
    if (flg!=PETSC_TRUE) {
      sprintf(str,"ilu");
      PetscOptionsSetValue(ostr,str);
      PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
    }

    /* by default use Manteuffel shift */
/*
#undef ostr
#define ostr "-psolve_pc_factor_shift_positive_definite"
    ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
    if (flg!=PETSC_TRUE) {
      sprintf(str,"1");
      PetscOptionsSetValue(ostr,str);
      PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
    }
*/
  }

  /* use prefix for options */
  ierr = KSPSetOptionsPrefix(ksp_Jac,"psolve_");CHKERRQ(ierr);

  /* possibly override defaults with settings from user defined options */
  ierr = KSPSetFromOptions(ksp_Jac);CHKERRQ(ierr);

/*
  ierr = KSPSetInitialGuessKnoll(ksp_Jac,PETSC_TRUE);CHKERRQ(ierr);
  ierr = KSPSetInitialGuessNonzero(ksp_Jac,PETSC_TRUE);CHKERRQ(ierr);
*/

  /* KSPSetUp only necessary to get more useful output from KSPView */
  /* FIXME: bug(?) in PETSc 2.2.1: for size==1 leads to error "Zero pivot row 0" */
  if (size>1) {
    ierr = KSPSetUp(ksp_Jac);CHKERRQ(ierr);
  }

  PetscPrintf(PETSC_COMM_WORLD,"KSP for preconditioning: P*z=r:\n");
#if defined(SUPERLU) || defined(BLOCKSOLVE)
  PetscReal rtol,atol,divtol;
  int maxits;
  ierr = KSPGetTolerances(ksp_Jac,&rtol,&atol,&divtol,&maxits);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"rtol: %g atol: %g divtol: %g maxits: %i\n",
    rtol,atol,divtol,maxits
  );
#else
  ierr = KSPView(ksp_Jac,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
#endif

  MagparFunctionLogReturn(0);
}


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)
{
  GridData     *gdata = (GridData*)P_data;
  MatStructure str = DIFFERENT_NONZERO_PATTERN;

  MagparFunctionInfoBegin;

  /* jok - TRUE means reuse current Jacobian else recompute Jacobian */
  if (jok) {
    ierr     = MatCopy(vode_pmat,vode_Jac,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
    str      = SAME_NONZERO_PATTERN;
    *jcurPtr = FALSE;
  } else {
    /* make PETSc vector tmpy point to PVODE vector y */
#ifdef UNIPROC
    ierr = VecPlaceArray(gdata->M,NV_DATA_S(y));CHKERRQ(ierr);
#else
    ierr = VecPlaceArray(gdata->M,NV_DATA_P(y));CHKERRQ(ierr);
#endif

    /* TODO: make Htot not calculate demag - we don't use it for Jacobian, could save (a little) time */
    ierr = Htot(gdata);CHKERRQ(ierr);

    /* compute the Jacobian */
    ierr = myLLGJacobian(PETSC_NULL,gdata->time,gdata->M,&vode_Jac,&vode_Jac,&str,gdata);CHKERRQ(ierr);

    ierr = VecResetArray(gdata->M);CHKERRQ(ierr);

    /* copy the Jacobian matrix */
    ierr = MatCopy(vode_Jac,vode_pmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);

    *jcurPtr = TRUE;
  }

  /* construct I-gamma*Jac  */

  ierr = MatScale(vode_Jac,-gamma);CHKERRQ(ierr);
  ierr = MatShift(vode_Jac,1.0);CHKERRQ(ierr);

  ierr = KSPSetOperators(ksp_Jac,vode_Jac,vode_Jac,str);CHKERRQ(ierr);

  MagparFunctionProfReturn(0);
}


/* adapted from $PETSC_DIR/src/ts/impls/implicit/pvode/petscpvode.c */

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)
{
  MagparFunctionInfoBegin;

  /* as defined in call to CVSpgmr in pvode_init: 1==left, 2==right */
  assert(lr==1);

#ifdef UNIPROC
  ierr = VecPlaceArray(vode_w1,NV_DATA_S(r));CHKERRQ(ierr);
  ierr = VecPlaceArray(vode_w2,NV_DATA_S(z));CHKERRQ(ierr);
#else
  ierr = VecPlaceArray(vode_w1,NV_DATA_P(r));CHKERRQ(ierr);
  ierr = VecPlaceArray(vode_w2,NV_DATA_P(z));CHKERRQ(ierr);
#endif

/* no difference in performance visible with mumag std prob. #4
   between abstol=0.1*delta*abstol and reltol=0.01
   because just 1 interation is usually sufficient
   (except if the Jacobian has been recalculated)

  KSP ksp_Jac;
  ierr = KSPSetTolerances(ksp_Jac,0.01,0.0,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);

  PetscInfo3(0,"set atol to: 0.1*delta*ewt[0]=%g*%g=%g\n",
    delta,NV_DATA_P(ewt)[0],0.1*delta*NV_DATA_P(ewt)[0]
  );
  ierr = KSPSetTolerances(ksp_Jac,0.0,0.1*delta*NV_DATA_P(ewt)[0],PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
*/

/* TODO: useful/necessary before every call to KSPSolve?
*/
/*
  ierr = KSPSetInitialGuessKnoll(ksp_Jac,PETSC_TRUE);CHKERRQ(ierr);
  ierr = KSPSetInitialGuessNonzero(ksp_Jac,PETSC_TRUE);CHKERRQ(ierr);
*/
  int its1;
  ierr = KSPSolve(ksp_Jac,vode_w1,vode_w2);CHKERRQ(ierr);

  ierr = VecResetArray(vode_w1);CHKERRQ(ierr);
  ierr = VecResetArray(vode_w2);CHKERRQ(ierr);

  ierr = KSPGetIterationNumber(ksp_Jac,(PetscInt*)&its1);CHKERRQ(ierr);
  KSPConvergedReason reason;
  ierr = KSPGetConvergedReason(ksp_Jac,&reason);CHKERRQ(ierr);

  ierr = PetscGetTime(&t_t2);CHKERRQ(ierr);
  PetscInfo2(0,"timing: %s took %g s\n",__FUNCT__,t_t2-t_t1);
  PetscInfo2(0,"psolve: %i its, reason: %i\n",its1,reason);
  if (reason<0) {
    PetscReal rnorm;

    ierr = KSPGetResidualNorm(ksp_Jac,&rnorm);CHKERRQ(ierr);
    PetscPrintf(PETSC_COMM_WORLD,
      "Warning: psolve: KSP did not converge - reason: %i, rnorm: %g, its: %i\n",
      reason,rnorm,its1
    );
  }

  MagparFunctionProfReturn(0);
}


int Jtimes(N_Vector v,N_Vector Jv,realtype t,N_Vector y,N_Vector fy,void *jac_data,N_Vector tmp)
{
  MagparFunctionInfoBegin;

#ifdef UNIPROC
  ierr = VecPlaceArray(vode_w1,NV_DATA_S(v));CHKERRQ(ierr);
  ierr = VecPlaceArray(vode_w2,NV_DATA_S(Jv));CHKERRQ(ierr);
#else
  ierr = VecPlaceArray(vode_w1,NV_DATA_P(v));CHKERRQ(ierr);
  ierr = VecPlaceArray(vode_w2,NV_DATA_P(Jv));CHKERRQ(ierr);
#endif

  ierr = MatMult(vode_pmat,vode_w1,vode_w2);CHKERRQ(ierr);

  ierr = VecResetArray(vode_w1);CHKERRQ(ierr);
  ierr = VecResetArray(vode_w2);CHKERRQ(ierr);

  MagparFunctionProfReturn(0);
}