PetscComputePreconMatImpl.cpp

Go to the documentation of this file.

#ifdef with_petsc

#include <stdio.h>
#include <arrayfunctions.h>
#include <simulation_object.h>
#include <mpivars_cpp.h>
#include <petscinterface.h>

#undef __FUNCT__
#define __FUNCT__ "PetscComputePreconMatImpl"

int PetscComputePreconMatImpl(  Mat Pmat,   
                                Vec Y,      
                                void *ctxt   )
{
  PetscErrorCode ierr;
  PETScContext* context = (PETScContext*) ctxt;
  SimulationObject* sim = (SimulationObject*) context->simobj;
  int nsims = context->nsims;

  PetscFunctionBegin;
  /* initialize preconditioning matrix to zero */
  MatZeroEntries(Pmat);

  /* copy solution from PETSc vector */
  for (int ns = 0; ns < nsims; ns++) {

    TransferVecFromPETSc( sim[ns].solver.u,
                          Y,
                          context,
                          ns,
                          context->offsets[ns]);

    HyPar* solver( &(sim[ns].solver) );
    MPIVariables* mpi( &(sim[ns].mpi) );

    int ndims = solver->ndims,
        nvars = solver->nvars,
        npoints = solver->npoints_local,
        ghosts = solver->ghosts,
        *dim = solver->dim_local,
        *isPeriodic = solver->isPeriodic,
        *points = context->points[ns],
        index[ndims],indexL[ndims],indexR[ndims],
        rows[nvars],cols[nvars];
  
    double *u = solver->u, 
           *iblank = solver->iblank,
           dxinv, values[nvars*nvars];
  
    /* apply boundary conditions and exchange data over MPI interfaces */
    solver->ApplyBoundaryConditions(solver,mpi,u,NULL,context->waqt);
    MPIExchangeBoundariesnD(ndims,nvars,dim,ghosts,mpi,u);
  
    /* loop through all computational points */
    for (int n = 0; n < npoints; n++) {
      int *this_point = points + n*(ndims+1);
      int p = this_point[ndims];
      int index[ndims]; _ArrayCopy1D_(this_point,index,ndims);
  
      double iblank = solver->iblank[p];
  
      /* compute the contributions from the hyperbolic flux derivatives along each dimension */
      if (solver->JFunction) {
        for (int dir = 0; dir < ndims; dir++) {
    
          /* compute indices and global 1D indices for left and right neighbors */
          _ArrayCopy1D_(index,indexL,ndims); indexL[dir]--;
          int pL;  _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
    
          _ArrayCopy1D_(index,indexR,ndims); indexR[dir]++;
          int pR;  _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
    
          int pg, pgL, pgR;
          pg  = (int) context->globalDOF[ns][p];
          pgL = (int) context->globalDOF[ns][pL];
          pgR = (int) context->globalDOF[ns][pR];
    
          /* Retrieve 1/delta-x at this grid point */
          _GetCoordinate_(dir,index[dir],dim,ghosts,solver->dxinv,dxinv);
    
          /* diagonal element */
          for (int v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pg + v; }
          solver->JFunction(values,(u+nvars*p),solver->physics,dir,nvars,0);
          _ArrayScale1D_(values,(dxinv*iblank),(nvars*nvars));
          MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES);
    
          /* left neighbor */
          if (pgL >= 0) {
            for (int v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pgL + v; }
            solver->JFunction(values,(u+nvars*pL),solver->physics,dir,nvars,1);
            _ArrayScale1D_(values,(-dxinv*iblank),(nvars*nvars));
            MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES);
          }
          
          /* right neighbor */
          if (pgR >= 0) {
            for (int v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pgR + v; }
            solver->JFunction(values,(u+nvars*pR),solver->physics,dir,nvars,-1);
            _ArrayScale1D_(values,(-dxinv*iblank),(nvars*nvars));
            MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES);
          }
        }
      }

      /* compute the contributions from the parabolic term derivatives along each dimension */
      if (solver->KFunction) {
        for (int dir = 0; dir < ndims; dir++) {
    
          /* compute indices and global 1D indices for left and right neighbors */
          _ArrayCopy1D_(index,indexL,ndims); indexL[dir]--;
          int pL;  _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
    
          _ArrayCopy1D_(index,indexR,ndims); indexR[dir]++;
          int pR;  _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
    
          int pg, pgL, pgR;
          pg  = (int) context->globalDOF[ns][p];
          pgL = (int) context->globalDOF[ns][pL];
          pgR = (int) context->globalDOF[ns][pR];
    
          /* Retrieve 1/delta-x at this grid point */
          _GetCoordinate_(dir,index[dir],dim,ghosts,solver->dxinv,dxinv);
   
          /* diagonal element */
          for (int v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pg + v; }
          solver->KFunction(values,(u+nvars*p),solver->physics,dir,nvars);
          _ArrayScale1D_(values,(-2*dxinv*dxinv*iblank),(nvars*nvars));
          MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES);
    
          /* left neighbor */
          if (pgL >= 0) {
            for (int v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pgL + v; }
            solver->KFunction(values,(u+nvars*pL),solver->physics,dir,nvars);
            _ArrayScale1D_(values,(dxinv*dxinv*iblank),(nvars*nvars));
            MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES);
          }
          
          /* right neighbor */
          if (pgR >= 0) {
            for (int v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pgR + v; }
            solver->KFunction(values,(u+nvars*pR),solver->physics,dir,nvars);
            _ArrayScale1D_(values,(dxinv*dxinv*iblank),(nvars*nvars));
            MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES);
          }
        }
      }
    }
  }

  MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd  (Pmat,MAT_FINAL_ASSEMBLY);
  
  MatShift(Pmat,context->shift);
  PetscFunctionReturn(0);
}

#endif