11 #ifndef _PETSCSOLVERI_H_ 12 #define _PETSCSOLVERI_H_ 19 #include "NamespaceHeader.H" 27 :m_homogeneous(false),
35 m_nz_init_guess(false),
110 const T& rhs = a_rhs;
121 const T& rhs = a_rhs;
133 CH_TIME(
"PetscSolver::create_mat_vec");
138 const PetscInt nc = a_phi.nComp();
140 MPI_Comm wcomm = Chombo_MPI::comm;
142 MPI_Comm wcomm = PETSC_COMM_SELF;
150 IntVect idghosts = a_phi.ghostVect();
153 MayDay::Error(
"PetscSolver<T>::create_mat_vec: No ghost cells in input LevelData<>.");
162 for ( dit = a_phi.dataIterator() ; dit.
ok() ; ++dit )
169 const PetscInt NN = nc*data;
173 PetscDataTypeToMPIDataType(PETSC_INT,&mtype);
174 MPI_Scan( &data, &result, 1, mtype, MPI_SUM, wcomm );
180 for ( dit = a_phi.dataIterator() ; dit.
ok() ; ++dit )
195 PetscInt *d_nnz=PETSC_NULL, *o_nnz=PETSC_NULL;
202 MPI_Comm wcomm = Chombo_MPI::comm;
204 PetscDataTypeToMPIDataType(PETSC_INT,&mtype);
205 MPI_Allreduce(&nn,&nglobal,1,mtype,MPI_SUM,wcomm);
209 ierr = PetscMalloc( (NN+1)*
sizeof(PetscInt), &d_nnz ); CHKERRQ(ierr);
210 ierr = PetscMalloc( (NN+1)*
sizeof(PetscInt), &o_nnz ); CHKERRQ(ierr);
211 for (PetscInt kk=0;kk<NN;kk++) d_nnz[kk] = o_nnz[kk] = 0;
216 for (PetscInt kk=0;kk<NN;kk++)
218 if (d_nnz[kk] > NN) d_nnz[kk] = NN;
219 if (o_nnz[kk] > nglobal-NN) o_nnz[kk] = nglobal-NN;
224 ierr = MatCreate(wcomm,&
m_mat);CHKERRQ(ierr);
228 ierr = MatSetSizes(
m_mat,NN,NN,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
229 ierr = MatSetBlockSize(
m_mat,nc);CHKERRQ(ierr);
230 ierr = MatSetType(
m_mat,MATAIJ);CHKERRQ(ierr);
232 ierr = MatSetFromOptions(
m_mat ); CHKERRQ(ierr);
233 ierr = MatSeqAIJSetPreallocation(
m_mat,nnzrow, d_nnz);CHKERRQ(ierr);
234 ierr = MatMPIAIJSetPreallocation(
m_mat,nnzrow, d_nnz, nnzrow/2, o_nnz);CHKERRQ(ierr);
235 ierr = MatSetOption(
m_mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE) ;CHKERRQ(ierr);
237 #if defined(PETSC_HAVE_HYPRE) 238 ierr = MatHYPRESetPreallocation(
m_mat,nnzrow, d_nnz, nnzrow/2, o_nnz);CHKERRQ(ierr);
243 ierr = PetscFree( d_nnz ); CHKERRQ(ierr);
244 ierr = PetscFree( o_nnz ); CHKERRQ(ierr);
252 ierr = VecDuplicate(
m_bb, &
m_rr ); CHKERRQ(ierr);
275 MPI_Comm wcomm = Chombo_MPI::comm;
277 MPI_Comm wcomm = PETSC_COMM_SELF;
295 #if PETSC_VERSION_GE(3,6,0) 296 strcat (str,
"pc_gamg_square_graph 20");
298 strcat (str,
"pc_gamg_square_graph true");
300 #if PETSC_VERSION_GE(3,7,0) 301 ierr = PetscOptionsInsertString(PETSC_NULL,str);CHKERRQ(ierr);
303 ierr = PetscOptionsInsertString(str);CHKERRQ(ierr);
308 PetscBool ism = PETSC_FALSE;
309 #if PETSC_VERSION_GE(3,7,0) 310 PetscOptionsGetBool(PETSC_NULL,
m_prestring,
"-ksp_monitor",&ism,PETSC_NULL);
312 PetscOptionsGetBool(
m_prestring,
"-ksp_monitor",&ism,PETSC_NULL);
316 ierr = SNESCreate( wcomm, &
m_snes ); CHKERRQ(ierr);
321 ierr = SNESSetFromOptions(
m_snes ); CHKERRQ(ierr);
322 ierr = SNESGetKSP(
m_snes, &ksp ); CHKERRQ(ierr);
325 ierr = KSPMonitorSet(ksp,
ksp_monitor_pout,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
331 KSPCreate( wcomm, &
m_ksp );
336 ierr = KSPSetFromOptions(
m_ksp);CHKERRQ(ierr);
341 #if PETSC_VERSION_GE(3,5,0) 344 ierr = KSPSetOperators(
m_ksp,
m_mat,
m_mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
346 ierr = KSPSetInitialGuessNonzero(
m_ksp,
m_nz_init_guess ? PETSC_TRUE : PETSC_FALSE );CHKERRQ(ierr);
353 PetscInt sz,ind,bs,n,m;
354 #if PETSC_VERSION_LT(3,4,0) & PETSC_VERSION_RELEASE 359 ierr = KSPGetPC( ksp, &pc ); CHKERRQ(ierr);
360 ierr = PCGetType( pc, &type ); CHKERRQ(ierr);
361 ierr = MatGetBlockSize(
m_mat, &bs ); CHKERRQ( ierr );
362 if ( strcmp(type,PCGAMG) == 0 && bs > 1 )
367 ierr = MatGetLocalSize(
m_mat, &m, &n ); CHKERRQ(ierr);
369 ierr = PetscMalloc( (sz+1)*
sizeof(PetscReal), &coords ); CHKERRQ(ierr);
370 for ( dit = a_phi.dataIterator(), ind = 0 ; dit.
ok() ; ++dit )
372 const Box &
box = a_phi.getBoxes()[dit];
377 for (PetscInt k=0; k<
CH_SPACEDIM; k++) coords[ind++] = (PetscReal)iv[k];
382 ierr = PetscFree( coords ); CHKERRQ(ierr);
394 ierr = MatZeroEntries(
m_mat ); CHKERRQ(ierr);
402 ierr = SNESGetKSP(
m_snes, &ksp ); CHKERRQ(ierr);
404 #if PETSC_VERSION_GE(3,5,0) 405 ierr = KSPSetOperators(ksp,
m_mat,
m_mat); CHKERRQ(ierr);
407 ierr = KSPSetOperators(ksp,
m_mat,
m_mat,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
427 const PetscInt nc = a_rhs.nComp();
429 MPI_Comm wcomm = Chombo_MPI::comm;
431 MPI_Comm wcomm = PETSC_COMM_SELF;
436 MPI_Barrier(Chombo_MPI::comm);
446 ierr = VecSetOption(
m_bb,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
447 ierr = VecSetOption(
m_xx,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
449 ierr = VecSet(
m_xx, 0.);CHKERRQ(ierr);
450 ierr = VecSet(
m_bb, 0.);CHKERRQ(ierr);
454 #pragma omp parallel for private(addV,ierr) 455 for (
int mybox=0;mybox<nbox; mybox++)
466 PetscInt mm, ki = nc*gidsFab(iv,0);
467 for (mm=0;mm<nc;mm++,ki++)
469 PetscScalar v = rhsFab(iv,mm);
479 ierr = VecSetValues(
m_bb,1,&ki,&v,INSERT_VALUES);
481 ierr = VecSetValues(
m_xx,1,&ki,&v,INSERT_VALUES);
486 ierr = VecAssemblyBegin(
m_bb );CHKERRQ(ierr);
487 ierr = VecAssemblyEnd(
m_bb );CHKERRQ(ierr);
490 ierr = VecAssemblyBegin(
m_xx );CHKERRQ(ierr);
491 ierr = VecAssemblyEnd(
m_xx );CHKERRQ(ierr);
498 ierr = MatNullSpaceCreate(wcomm,PETSC_TRUE,0,PETSC_NULL,&nullsp);CHKERRQ(ierr);
500 ierr = MatSetNullSpace(
m_mat, nullsp );CHKERRQ(ierr);
505 MPI_Barrier(Chombo_MPI::comm);
531 CH_TIME(
"PetscSolver::apply_mfree");
535 MatShellGetContext(A, &ctx);
540 PetscFunctionReturn(0);
548 CH_TIME(
"PetscSolver::solve_mfree_private");
551 MPI_Comm wcomm = Chombo_MPI::comm;
553 MPI_Comm wcomm = PETSC_COMM_SELF;
561 ierr = MatGetSize(
m_mat, &M, &N);CHKERRQ(ierr);
CH_assert(M == N);
562 ierr = MatGetLocalSize(
m_mat, &m, &n);CHKERRQ(ierr);
564 ierr = MatCreateShell(wcomm,m,n,N,N,(
void *)
this,&L);CHKERRQ(ierr);
565 ierr = MatShellSetOperation(L,MATOP_MULT,(
void(*)(
void))
apply_mfree);
576 ierr = VecSetOption(
m_bb,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
577 ierr = VecSetOption(
m_xx,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
578 ierr = VecSet(
m_xx, 0.);CHKERRQ(ierr);
579 ierr = VecSet(
m_bb, 0.);CHKERRQ(ierr);
581 const PetscInt nc = a_rhs.nComp();
585 for (
int mybox=0;mybox<nbox; mybox++)
596 PetscInt mm, ki = nc*gidsFab(iv,0);
597 for (mm=0;mm<nc;mm++,ki++)
599 PetscScalar v = rhsFab(iv,mm);
609 ierr = VecSetValues(
m_bb,1,&ki,&v,INSERT_VALUES); CHKERRQ(ierr);
611 ierr = VecSetValues(
m_xx,1,&ki,&v,INSERT_VALUES); CHKERRQ(ierr);
616 ierr = VecAssemblyBegin(
m_bb );CHKERRQ(ierr);
617 ierr = VecAssemblyEnd(
m_bb );CHKERRQ(ierr);
618 ierr = VecAssemblyBegin(
m_xx );CHKERRQ(ierr);
619 ierr = VecAssemblyEnd(
m_xx );CHKERRQ(ierr);
622 KSPCreate( wcomm, &ksp );
627 ierr = KSPSetFromOptions( ksp );CHKERRQ(ierr);
629 #if PETSC_VERSION_GE(3,7,0) 630 PetscOptionsGetBool(PETSC_NULL,
m_prestring,
"-ksp_monitor",&ism,PETSC_NULL);
632 PetscOptionsGetBool(
m_prestring,
"-ksp_monitor",&ism,PETSC_NULL);
637 ierr = KSPMonitorSet(ksp,
ksp_monitor_pout,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
639 #if PETSC_VERSION_GE(3,5,0) 640 ierr = KSPSetOperators(ksp,L,
m_mat); CHKERRQ(ierr);
642 ierr = KSPSetOperators(ksp,L,
m_mat,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
645 ierr = KSPSetInitialGuessNonzero(ksp,
m_nz_init_guess ? PETSC_TRUE : PETSC_FALSE ); CHKERRQ(ierr);
648 PC pc; PetscInt sz,ind,bs,n,m;
649 #if PETSC_VERSION_LT(3,4,0) & PETSC_VERSION_RELEASE 654 ierr = KSPGetPC( ksp, &pc ); CHKERRQ(ierr);
655 ierr = PCGetType( pc, &type ); CHKERRQ(ierr);
656 ierr = MatGetBlockSize(
m_mat, &bs ); CHKERRQ( ierr );
657 if ( strcmp(type,PCGAMG) == 0 && bs > 1 )
662 ierr = MatGetLocalSize(
m_mat, &m, &n ); CHKERRQ(ierr);
664 ierr = PetscMalloc( (sz+1)*
sizeof(PetscReal), &coords ); CHKERRQ(ierr);
665 for ( dit = a_phi.dataIterator(), ind = 0 ; dit.
ok() ; ++dit )
667 const Box &
box = a_phi.getBoxes()[dit];
672 for (PetscInt k=0; k<
CH_SPACEDIM; k++) coords[ind++] = (PetscReal)iv[k];
677 ierr = PetscFree( coords ); CHKERRQ(ierr);
682 ierr = KSPSolve( ksp,
m_bb,
m_xx );CHKERRQ(ierr);
702 const PetscScalar *arr;
704 const PetscInt nc = a_phi.nComp();
706 ierr = VecGetArrayRead(xx,&arr); CHKERRQ(ierr);
710 for (
int mybox=0;mybox<nbox; mybox++)
720 PetscInt mm, ki = nc*gidsFab(iv,0);
721 for (mm=0;mm<nc;mm++,ki++)
723 phiFab(iv,mm) = arr[ki - nc*
m_gid0];
727 ierr = VecRestoreArrayRead(xx,&arr); CHKERRQ(ierr);
735 CH_TIME(
"PetscSolver::putChomboInPetsc");
739 const PetscInt nc = a_phi.nComp();
743 ierr = VecSet( out, 0.); CHKERRQ(ierr);
747 for (
int mybox=0;mybox<nbox; mybox++)
757 PetscInt mm, ki = nc*gidsFab(iv,0);
758 for (mm=0;mm<nc;mm++,ki++)
760 PetscScalar v = phiFab(iv,mm);
762 if (std::isnan(v)) v = 0;
763 ierr = VecSetValues(out,1,&ki,&v,INSERT_VALUES);
768 ierr = VecAssemblyBegin( out ); CHKERRQ(ierr);
769 ierr = VecAssemblyEnd( out ); CHKERRQ(ierr);
783 PetscScalar alpha = -1;
786 ierr = VecDuplicate(
m_xx, &tempVec);CHKERRQ(ierr);
787 ierr = MatMult(
m_mat,
m_xx, tempVec);CHKERRQ(ierr);
788 ierr = VecAXPY(tempVec,alpha,
m_bb);CHKERRQ(ierr);
792 Residual.
define( dbl, 0, idghosts );
795 VecDestroy( &tempVec );
813 ierr = VecNorm(
m_rr,NORM_INFINITY,&norm);CHKERRQ(ierr);
844 #include "NamespaceFooter.H" bool ok()
Definition: BoxIterator.H:281
#define CH_TIMERS(name)
Definition: CH_Timer.H:133
Definition: LinearSolver.H:28
virtual void define(Real a_dx, bool a_homogeneous=true)
Definition: PetscSolverI.H:81
#define CH_SPACEDIM
Definition: SPACE.H:51
#define CH_assert(cond)
Definition: CHArray.H:37
KSP m_ksp
Definition: PetscSolver.H:243
virtual void solve(T &a_phi, const T &a_rhs)
Definition: PetscSolverI.H:106
int solve_private(T &a_phi, const T &a_rhs)
Definition: PetscSolverI.H:417
PetscInt m_gid0
Definition: PetscSolver.H:263
void setVal(T a_x, const Box &a_bx, int a_nstart, int a_numcomp)
{ data modification functions}
Definition: BaseFabImplem.H:399
Real norm(const BoxLayoutData< FArrayBox > &A, const Interval &interval, const int &p=2)
#define CH_START(tpointer)
Definition: CH_Timer.H:145
static PetscErrorCode ksp_monitor_pout(KSP ksp, PetscInt it, PetscReal rnorm, void *ctx)
Definition: PetscSolver.H:65
virtual int formMatrix(Mat, const T *=0, PetscInt my0=0, PetscInt nloc=0, PetscInt *d=0, PetscInt *o=0)=0
virtual bool ok() const
return true if this iterator is still in its Layout
Definition: LayoutIterator.H:117
Definition: DataIterator.H:190
virtual void solve_mfree(T &a_phi, const T &a_rhs, LinearOp< T > *a_op)
Definition: PetscSolverI.H:116
iterates through the IntVects of a Box
Definition: BoxIterator.H:37
virtual Real addBCrhsValue(const IntVect &a_iv, const T &a_phi, const DataIndex &a_datInd, const Real &coeff=1)
Definition: PetscSolver.H:230
virtual void exchange(void)
Simplest case – do all components.
Definition: LevelDataI.H:451
#define CH_TIMER(name, tpointer)
Definition: CH_Timer.H:63
Definition: EBInterface.H:45
virtual int applyOp(T &a_phi, const T &a_rhs)
Definition: PetscSolverI.H:823
virtual ~PetscSolver()
Definition: PetscSolverI.H:74
virtual BaseFab< Real > & getRegFab(T &a_fab, const DataIndex &a_datInd)=0
Real m_dx
Definition: PetscSolver.H:194
virtual bool supportNNZExact() const
Definition: PetscSolver.H:173
bool m_null
Definition: PetscSolver.H:255
void const char const int const int const int const Real const Real * A
Definition: Lapack.H:83
int size() const
Definition: DataIterator.H:218
Real computeResidual()
Definition: PetscSolverI.H:779
PetscErrorCode putPetscInChombo(T &a_phi, const Vec xx)
Definition: PetscSolverI.H:699
void setNull(bool n=true)
Definition: PetscSolverI.H:99
LevelData< BaseFab< bool > > m_bccode
Definition: PetscSolver.H:257
LinearOp< T > * m_op_mfree
Definition: PetscSolver.H:110
#define CH_TIME(name)
Definition: CH_Timer.H:82
int solve_mfree_private(T &a_phi, const T &a_rhs, LinearOp< T > *a_op)
Definition: PetscSolverI.H:544
Mat m_mat
Definition: PetscSolver.H:238
Vec m_bb
Definition: PetscSolver.H:241
PetscErrorCode(* m_function)(SNES, Vec, Vec, void *)
Definition: PetscSolver.H:246
double Real
Definition: REAL.H:33
virtual void define(const DisjointBoxLayout &dp, int comps, const IntVect &ghost=IntVect::Zero, const DataFactory< T > &a_factory=DefaultDataFactory< T >())
Definition: LevelDataI.H:80
LevelData< BaseFab< PetscInt > > m_gids
Definition: PetscSolver.H:262
A BoxLayout that has a concept of disjointedness.
Definition: DisjointBoxLayout.H:30
Box get(const LayoutIndex &it) const
Definition: BoxLayout.H:717
const IntVect & ghostVect() const
Definition: LevelData.H:170
static void Error(const char *const a_msg=m_nullString, int m_exitCode=CH_DEFAULT_ERROR_CODE)
Print out message to cerr and exit with the specified exit code.
virtual int getNNZPerRow() const
Definition: PetscSolver.H:166
Real normInfinity(const T &a_phi) const
Definition: PetscSolverI.H:807
PetscSolver()
Definition: PetscSolverI.H:26
static const IntVect Zero
Definition: IntVect.H:658
virtual void clear(T &a_lhs)
Definition: LinearSolver.H:70
T m_phi_mfree
Definition: PetscSolver.H:111
bool m_nz_init_guess
Definition: PetscSolver.H:256
char m_prestring[32]
Definition: PetscSolver.H:254
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:469
static PetscErrorCode apply_mfree(Mat A, Vec x, Vec f)
Definition: PetscSolverI.H:529
virtual void create(T &a_lhs, const T &a_rhs)=0
Vec m_xx
Definition: PetscSolver.H:241
Definition: DataIndex.H:114
T m_Lphi_mfree
Definition: PetscSolver.H:112
const DisjointBoxLayout & disjointBoxLayout() const
Definition: LevelData.H:209
void begin()
Definition: BoxIterator.H:150
#define CH_STOP(tpointer)
Definition: CH_Timer.H:150
int setup_solver(const T &a_phi)
Definition: PetscSolverI.H:264
PetscErrorCode putChomboInPetsc(Vec out, const T &a_phi)
Definition: PetscSolverI.H:733
Definition: PetscSolver.H:36
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
void destroy()
Definition: PetscSolverI.H:43
PetscErrorCode(* m_jacobian)(SNES, Vec, Mat *, Mat *, MatStructure *, void *)
Definition: PetscSolver.H:250
void const char const int * M
Definition: Lapack.H:83
void next()
Definition: BoxIterator.H:168
void const char const int const int * N
Definition: Lapack.H:83
SNES m_snes
Definition: PetscSolver.H:242
int create_mat_vec(const T &a_phi)
Definition: PetscSolverI.H:131
PetscInt m_defined
Definition: PetscSolver.H:245
virtual Real dx() const
Definition: LinearSolver.H:128
Vec m_rr
Definition: PetscSolver.H:241
bool m_homogeneous
Definition: PetscSolver.H:188