本文整理汇总了C++中KSPGetPC函数的典型用法代码示例。如果您正苦于以下问题：C++ KSPGetPC函数的具体用法？C++ KSPGetPC怎么用？C++ KSPGetPC使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了KSPGetPC函数的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: PetscConvEstSetSolver
//.........这里部分代码省略.........
ierr = DMSetCoarseDM(dm[r], dm[r-1]);CHKERRQ(ierr);
ierr = DMCopyDisc(ce->idm, dm[r]);CHKERRQ(ierr);
ierr = DMCopyTransform(ce->idm, dm[r]);CHKERRQ(ierr);
ierr = PetscObjectGetName((PetscObject) dm[r-1], &dmname);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) dm[r], dmname);CHKERRQ(ierr);
for (f = 0; f <= ce->Nf; ++f) {
PetscErrorCode (*nspconstr)(DM, PetscInt, MatNullSpace *);
ierr = DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);CHKERRQ(ierr);
ierr = DMSetNullSpaceConstructor(dm[r], f, nspconstr);CHKERRQ(ierr);
}
}
ierr = DMViewFromOptions(dm[r], NULL, "-conv_dm_view");CHKERRQ(ierr);
/* Create solution */
ierr = DMCreateGlobalVector(dm[r], &u);CHKERRQ(ierr);
ierr = DMGetField(dm[r], 0, NULL, &disc);CHKERRQ(ierr);
ierr = PetscObjectGetName(disc, &uname);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) u, uname);CHKERRQ(ierr);
/* Setup solver */
ierr = SNESReset(ce->snes);CHKERRQ(ierr);
ierr = SNESSetDM(ce->snes, dm[r]);CHKERRQ(ierr);
ierr = DMPlexSetSNESLocalFEM(dm[r], ctx, ctx, ctx);CHKERRQ(ierr);
ierr = SNESSetFromOptions(ce->snes);CHKERRQ(ierr);
/* Create initial guess */
ierr = DMProjectFunction(dm[r], t, ce->initGuess, ce->ctxs, INSERT_VALUES, u);CHKERRQ(ierr);
ierr = SNESSolve(ce->snes, NULL, u);CHKERRQ(ierr);
ierr = PetscLogEventBegin(event, ce, 0, 0, 0);CHKERRQ(ierr);
ierr = DMComputeL2FieldDiff(dm[r], t, ce->exactSol, ce->ctxs, u, &ce->errors[r*ce->Nf]);CHKERRQ(ierr);
ierr = PetscLogEventEnd(event, ce, 0, 0, 0);CHKERRQ(ierr);
for (f = 0; f < ce->Nf; ++f) {
PetscSection s, fs;
PetscInt lsize;
/* Could use DMGetOutputDM() to add in Dirichlet dofs */
ierr = DMGetSection(dm[r], &s);CHKERRQ(ierr);
ierr = PetscSectionGetField(s, f, &fs);CHKERRQ(ierr);
ierr = PetscSectionGetConstrainedStorageSize(fs, &lsize);CHKERRQ(ierr);
ierr = MPI_Allreduce(&lsize, &dof[r*ce->Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) ce->snes));CHKERRQ(ierr);
ierr = PetscLogEventSetDof(event, f, dof[r*ce->Nf+f]);CHKERRQ(ierr);
ierr = PetscLogEventSetError(event, f, ce->errors[r*ce->Nf+f]);CHKERRQ(ierr);
}
/* Monitor */
if (ce->monitor) {
PetscReal *errors = &ce->errors[r*ce->Nf];
ierr = PetscPrintf(comm, "L_2 Error: ");CHKERRQ(ierr);
if (ce->Nf > 1) {ierr = PetscPrintf(comm, "[");CHKERRQ(ierr);}
for (f = 0; f < ce->Nf; ++f) {
if (f > 0) {ierr = PetscPrintf(comm, ", ");CHKERRQ(ierr);}
if (errors[f] < 1.0e-11) {ierr = PetscPrintf(comm, "< 1e-11");CHKERRQ(ierr);}
else {ierr = PetscPrintf(comm, "%g", (double)errors[f]);CHKERRQ(ierr);}
}
if (ce->Nf > 1) {ierr = PetscPrintf(comm, "]");CHKERRQ(ierr);}
ierr = PetscPrintf(comm, "\n");CHKERRQ(ierr);
}
if (!r) {
/* PCReset() does not wipe out the level structure */
KSP ksp;
PC pc;
ierr = SNESGetKSP(ce->snes, &ksp);CHKERRQ(ierr);
ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
ierr = PCMGGetLevels(pc, &oldnlev);CHKERRQ(ierr);
}
/* Cleanup */
ierr = VecDestroy(&u);CHKERRQ(ierr);
ierr = PetscLogStagePop();CHKERRQ(ierr);
}
for (r = 1; r <= Nr; ++r) {
ierr = DMDestroy(&dm[r]);CHKERRQ(ierr);
}
/* Fit convergence rate */
ierr = PetscMalloc2(Nr+1, &x, Nr+1, &y);CHKERRQ(ierr);
for (f = 0; f < ce->Nf; ++f) {
for (r = 0; r <= Nr; ++r) {
x[r] = PetscLog10Real(dof[r*ce->Nf+f]);
y[r] = PetscLog10Real(ce->errors[r*ce->Nf+f]);
}
ierr = PetscLinearRegression(Nr+1, x, y, &slope, &intercept);CHKERRQ(ierr);
/* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
alpha[f] = -slope * dim;
}
ierr = PetscFree2(x, y);CHKERRQ(ierr);
ierr = PetscFree2(dm, dof);CHKERRQ(ierr);
/* Restore solver */
ierr = SNESReset(ce->snes);CHKERRQ(ierr);
{
/* PCReset() does not wipe out the level structure */
KSP ksp;
PC pc;
ierr = SNESGetKSP(ce->snes, &ksp);CHKERRQ(ierr);
ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
ierr = PCMGSetLevels(pc, oldnlev, NULL);CHKERRQ(ierr);
ierr = DMSetRefineLevel(ce->idm, oldlevel);CHKERRQ(ierr); /* The damn DMCoarsen() calls in PCMG can reset this */
}
ierr = SNESSetDM(ce->snes, ce->idm);CHKERRQ(ierr);
ierr = DMPlexSetSNESLocalFEM(ce->idm, ctx, ctx, ctx);CHKERRQ(ierr);
ierr = SNESSetFromOptions(ce->snes);CHKERRQ(ierr);
PetscFunctionReturn(0);
}示例2: main
//.........这里部分代码省略.........
ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
for (Ii=rstart; Ii<rend; Ii++) {
v = -1.0; i = Ii/n; j = Ii - i*n;
if (i>0) {
J = Ii-n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
if (i<n-1) {
J = Ii+n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
if (j>0) {
J = Ii-1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
if (j<n-1) {
J = Ii+1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
v = 4.0 - sigma1*h2;
ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* Check whether A is symmetric */
ierr = PetscOptionsHasName(PETSC_NULL, "-check_symmetric", &flg);CHKERRQ(ierr);
if (flg) {
Mat Trans;
ierr = MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);
ierr = MatEqual(A, Trans, &flg);
if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"A is not symmetric");
ierr = MatDestroy(&Trans);CHKERRQ(ierr);
}
ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
/* make A complex Hermitian */
Ii = 0; J = dim-1;
if (Ii >= rstart && Ii < rend){
v = sigma2*h2; /* RealPart(v) = 0.0 */
ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);
v = -sigma2*h2;
ierr = MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr);
}
Ii = dim-2; J = dim-1;
if (Ii >= rstart && Ii < rend){
v = sigma2*h2; /* RealPart(v) = 0.0 */
ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);
v = -sigma2*h2;
ierr = MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* Check whether A is Hermitian */
ierr = PetscOptionsHasName(PETSC_NULL, "-check_Hermitian", &flg);CHKERRQ(ierr);
if (flg) {
Mat Hermit;
if (disp_mat){
if (!rank) printf(" A:\n");
ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
}
ierr = MatHermitianTranspose(A,MAT_INITIAL_MATRIX, &Hermit);
if (disp_mat){
if (!rank) printf(" A_Hermitian:\n");
ierr = MatView(Hermit,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
}
ierr = MatEqual(A, Hermit, &flg);
if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"A is not Hermitian");
ierr = MatDestroy(&Hermit);CHKERRQ(ierr);
}
ierr = MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);CHKERRQ(ierr);
/* Create a Hermitian matrix As in sbaij format */
ierr = MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&As);CHKERRQ(ierr);
if (disp_mat){
if (!rank) {ierr = PetscPrintf(PETSC_COMM_SELF," As:\n");CHKERRQ(ierr);}
ierr = MatView(As,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
}
/* Test MatGetInertia() */
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,As,As,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCCHOLESKY);CHKERRQ(ierr);
ierr = PCSetFromOptions(pc);CHKERRQ(ierr);
ierr = PCSetUp(pc);CHKERRQ(ierr);
ierr = PCFactorGetMatrix(pc,&F);CHKERRQ(ierr);
ierr = MatGetInertia(F,&nneg,&nzero,&npos);CHKERRQ(ierr);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
if (!rank){
ierr = PetscPrintf(PETSC_COMM_SELF," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);CHKERRQ(ierr);
}
/* Free spaces */
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
if (use_random) {ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);}
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = MatDestroy(&As);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}示例3: main
int main(int argc,char **argv)
{
PetscErrorCode ierr;
PetscInt its,n,Nx=PETSC_DECIDE,Ny=PETSC_DECIDE,nlocal,i;
PetscMPIInt size;
PC pc;
PetscInt mx,my;
Mat A;
GridCtx fine_ctx;
KSP ksp;
PetscBool flg;
PetscInitialize(&argc,&argv,NULL,help);
/* set up discretization matrix for fine grid */
/* ML requires input of fine-grid matrix. It determines nlevels. */
fine_ctx.mx = 9; fine_ctx.my = 9;
ierr = PetscOptionsGetInt(NULL,"-mx",&mx,&flg);CHKERRQ(ierr);
if (flg) fine_ctx.mx = mx;
ierr = PetscOptionsGetInt(NULL,"-my",&my,&flg);CHKERRQ(ierr);
if (flg) fine_ctx.my = my;
ierr = PetscPrintf(PETSC_COMM_WORLD,"Fine grid size %D by %D\n",fine_ctx.mx,fine_ctx.my);CHKERRQ(ierr);
n = fine_ctx.mx*fine_ctx.my;
MPI_Comm_size(PETSC_COMM_WORLD,&size);
ierr = PetscOptionsGetInt(NULL,"-Nx",&Nx,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,"-Ny",&Ny,NULL);CHKERRQ(ierr);
ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,fine_ctx.mx,
fine_ctx.my,Nx,Ny,1,1,NULL,NULL,&fine_ctx.da);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(fine_ctx.da,&fine_ctx.x);CHKERRQ(ierr);
ierr = VecDuplicate(fine_ctx.x,&fine_ctx.b);CHKERRQ(ierr);
ierr = VecGetLocalSize(fine_ctx.x,&nlocal);CHKERRQ(ierr);
ierr = DMCreateLocalVector(fine_ctx.da,&fine_ctx.localX);CHKERRQ(ierr);
ierr = VecDuplicate(fine_ctx.localX,&fine_ctx.localF);CHKERRQ(ierr);
ierr = MatCreateAIJ(PETSC_COMM_WORLD,nlocal,nlocal,n,n,5,NULL,3,NULL,&A);CHKERRQ(ierr);
ierr = FormJacobian_Grid(&fine_ctx,&A);CHKERRQ(ierr);
/* create linear solver */
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCML);CHKERRQ(ierr);
/* set options, then solve system */
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); /* calls PCSetFromOptions_MG/ML */
for (i=0; i<3; i++) {
if (i<2) { /* test DIFFERENT_NONZERO_PATTERN */
/* set values for rhs vector */
ierr = VecSet(fine_ctx.b,i+1.0);CHKERRQ(ierr);
/* modify A */
ierr = MatShift(A,1.0);CHKERRQ(ierr);
ierr = MatScale(A,2.0);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
} else { /* test SAME_NONZERO_PATTERN */
ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
}
ierr = KSPSolve(ksp,fine_ctx.b,fine_ctx.x);CHKERRQ(ierr);
ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %D\n",its);CHKERRQ(ierr);
}
/* free data structures */
ierr = VecDestroy(&fine_ctx.x);CHKERRQ(ierr);
ierr = VecDestroy(&fine_ctx.b);CHKERRQ(ierr);
ierr = DMDestroy(&fine_ctx.da);CHKERRQ(ierr);
ierr = VecDestroy(&fine_ctx.localX);CHKERRQ(ierr);
ierr = VecDestroy(&fine_ctx.localF);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}示例4: TaoSolve_NTR
static PetscErrorCode TaoSolve_NTR(Tao tao)
{
TAO_NTR *tr = (TAO_NTR *)tao->data;
PC pc;
KSPConvergedReason ksp_reason;
TaoConvergedReason reason;
PetscReal fmin, ftrial, prered, actred, kappa, sigma, beta;
PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius;
PetscReal f, gnorm;
PetscReal delta;
PetscReal norm_d;
PetscErrorCode ierr;
PetscInt iter = 0;
PetscInt bfgsUpdates = 0;
PetscInt needH;
PetscInt i_max = 5;
PetscInt j_max = 1;
PetscInt i, j, N, n, its;
PetscFunctionBegin;
if (tao->XL || tao->XU || tao->ops->computebounds) {
ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by ntr algorithm\n");CHKERRQ(ierr);
}
tao->trust = tao->trust0;
/* Modify the radius if it is too large or small */
tao->trust = PetscMax(tao->trust, tr->min_radius);
tao->trust = PetscMin(tao->trust, tr->max_radius);
if (NTR_PC_BFGS == tr->pc_type && !tr->M) {
ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr);
ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr);
ierr = MatCreateLMVM(((PetscObject)tao)->comm,n,N,&tr->M);CHKERRQ(ierr);
ierr = MatLMVMAllocateVectors(tr->M,tao->solution);CHKERRQ(ierr);
}
/* Check convergence criteria */
ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr);
ierr = VecNorm(tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
needH = 1;
ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason);CHKERRQ(ierr);
if (reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
/* Create vectors for the limited memory preconditioner */
if ((NTR_PC_BFGS == tr->pc_type) &&
(BFGS_SCALE_BFGS != tr->bfgs_scale_type)) {
if (!tr->Diag) {
ierr = VecDuplicate(tao->solution, &tr->Diag);CHKERRQ(ierr);
}
}
switch(tr->ksp_type) {
case NTR_KSP_NASH:
ierr = KSPSetType(tao->ksp, KSPNASH);CHKERRQ(ierr);
if (tao->ksp->ops->setfromoptions) {
(*tao->ksp->ops->setfromoptions)(tao->ksp);
}
break;
case NTR_KSP_STCG:
ierr = KSPSetType(tao->ksp, KSPSTCG);CHKERRQ(ierr);
if (tao->ksp->ops->setfromoptions) {
(*tao->ksp->ops->setfromoptions)(tao->ksp);
}
break;
default:
ierr = KSPSetType(tao->ksp, KSPGLTR);CHKERRQ(ierr);
if (tao->ksp->ops->setfromoptions) {
(*tao->ksp->ops->setfromoptions)(tao->ksp);
}
break;
}
/* Modify the preconditioner to use the bfgs approximation */
ierr = KSPGetPC(tao->ksp, &pc);CHKERRQ(ierr);
switch(tr->pc_type) {
case NTR_PC_NONE:
ierr = PCSetType(pc, PCNONE);CHKERRQ(ierr);
if (pc->ops->setfromoptions) {
(*pc->ops->setfromoptions)(pc);
}
break;
case NTR_PC_AHESS:
ierr = PCSetType(pc, PCJACOBI);CHKERRQ(ierr);
if (pc->ops->setfromoptions) {
(*pc->ops->setfromoptions)(pc);
}
ierr = PCJacobiSetUseAbs(pc);CHKERRQ(ierr);
break;
case NTR_PC_BFGS:
ierr = PCSetType(pc, PCSHELL);CHKERRQ(ierr);
//.........这里部分代码省略.........
示例5: PETScMGSolver_UpdateSolvers
void PETScMGSolver_UpdateSolvers( PETScMGSolver* self ) {
PETScMGSolver_Level* level;
PC pc;
KSP levelKSP;
PC levelPC;
PetscErrorCode ec;
unsigned l_i;
PetscTruth smoothers_differ, flag;
PetscMPIInt size;
MPI_Comm comm;
assert( self && Stg_CheckType( self, PETScMGSolver ) );
ec = KSPGetPC( self->mgData->ksp, &pc );
CheckPETScError( ec );
ec = PCMGSetLevels( pc, self->nLevels, PETSC_NULL );
CheckPETScError( ec );
ec = PCMGSetType( pc, PC_MG_MULTIPLICATIVE );
CheckPETScError( ec );
ec=PetscOptionsGetTruth( PETSC_NULL, "-pc_mg_different_smoothers", &smoothers_differ, &flag ); CheckPETScError(ec);
ec=PetscObjectGetComm( (PetscObject)pc, &comm ); CheckPETScError(ec);
MPI_Comm_size( comm, &size );
for( l_i = 1; l_i < self->nLevels; l_i++ ) {
level = self->levels + l_i;
printf("Configuring MG level %d \n", l_i );
ec = PCMGGetSmootherDown( pc, l_i, &levelKSP );
CheckPETScError( ec );
if(smoothers_differ==PETSC_TRUE) { ec=KSPAppendOptionsPrefix( levelKSP, "down_" ); CheckPETScError(ec); }
ec = KSPSetType( levelKSP, KSPRICHARDSON ); CheckPETScError( ec );
ec = KSPGetPC( levelKSP, &levelPC ); CheckPETScError( ec );
if(size==1) {
ec = PCSetType( levelPC, PCSOR ); CheckPETScError( ec );
}
/* This does not work - bug with the order the operators are created I guess */
/* For parallel jobs you best bet is to use the command line args and let petsc work it out */
/*
else {
KSP *sub_ksp;
PetscInt k, n_local, first_local;
PC sub_pc;
PCSetType( levelPC, PCBJACOBI );
KSPSetUp( levelKSP );
PCBJacobiGetSubKSP( levelPC, &n_local,&first_local,&sub_ksp);
for(k=0;k<n_local;k++ ) {
KSPSetType( sub_ksp[k], KSPFGMRES );
KSPGetPC( sub_ksp[k], &sub_pc );
PCSetType( sub_pc, PCSOR );
}
}
*/
ec = KSPSetTolerances( levelKSP, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, level->nDownIts ); CheckPETScError( ec );
if( l_i == self->nLevels - 1 ) {
ec = KSPSetInitialGuessNonzero( levelKSP, PETSC_TRUE ); CheckPETScError( ec );
}
else { ec = KSPSetInitialGuessNonzero( levelKSP, PETSC_FALSE ); CheckPETScError( ec ); }
ec = PCMGGetSmootherUp( pc, l_i, &levelKSP ); CheckPETScError( ec );
if(smoothers_differ==PETSC_TRUE) { ec=KSPAppendOptionsPrefix( levelKSP, "up_" ); CheckPETScError(ec); }
ec = KSPSetType( levelKSP, KSPRICHARDSON ); CheckPETScError( ec );
ec = KSPGetPC( levelKSP, &levelPC ); CheckPETScError( ec );
if(size==1) {
ec = PCSetType( levelPC, PCSOR ); CheckPETScError( ec );
}
ec = KSPSetTolerances( levelKSP, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, level->nUpIts ); CheckPETScError( ec );
ec = KSPSetInitialGuessNonzero( levelKSP, PETSC_TRUE ); CheckPETScError( ec );
ec = PCMGSetCyclesOnLevel( pc, l_i, level->nCycles ); CheckPETScError( ec );
}
}示例6: main
//.........这里部分代码省略.........
idx2[7*i+4] = 9*i+6; idx2[7*i+5] = 9*i+7; idx2[7*i+6] = 9*i+8;
}
ierr = ISCreateGeneral(PETSC_COMM_WORLD,7*ngen,idx2,PETSC_COPY_VALUES,&user.is_diff);CHKERRQ(ierr);
ierr = ISComplement(user.is_diff,0,user.neqs_pgrid,&user.is_alg);CHKERRQ(ierr);
ierr = PetscFree(idx2);CHKERRQ(ierr);
/* Set run time options */
ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Transient stability fault options","");CHKERRQ(ierr);
{
user.tfaulton = 1.0;
user.tfaultoff = 1.2;
user.Rfault = 0.0001;
user.faultbus = 8;
ierr = PetscOptionsReal("-tfaulton","","",user.tfaulton,&user.tfaulton,NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-tfaultoff","","",user.tfaultoff,&user.tfaultoff,NULL);CHKERRQ(ierr);
ierr = PetscOptionsInt("-faultbus","","",user.faultbus,&user.faultbus,NULL);CHKERRQ(ierr);
user.t0 = 0.0;
user.tmax = 1.5;
ierr = PetscOptionsReal("-t0","","",user.t0,&user.t0,NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-tmax","","",user.tmax,&user.tmax,NULL);CHKERRQ(ierr);
user.freq_u = 61.0;
user.freq_l = 59.0;
user.pow = 2;
ierr = PetscOptionsReal("-frequ","","",user.freq_u,&user.freq_u,NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-freql","","",user.freq_l,&user.freq_l,NULL);CHKERRQ(ierr);
ierr = PetscOptionsInt("-pow","","",user.pow,&user.pow,NULL);CHKERRQ(ierr);
}
ierr = PetscOptionsEnd();CHKERRQ(ierr);
/* Create DMs for generator and network subsystems */
ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_gen,1,1,NULL,&user.dmgen);CHKERRQ(ierr);
ierr = DMSetOptionsPrefix(user.dmgen,"dmgen_");CHKERRQ(ierr);
ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_net,1,1,NULL,&user.dmnet);CHKERRQ(ierr);
ierr = DMSetOptionsPrefix(user.dmnet,"dmnet_");CHKERRQ(ierr);
/* Create a composite DM packer and add the two DMs */
ierr = DMCompositeCreate(PETSC_COMM_WORLD,&user.dmpgrid);CHKERRQ(ierr);
ierr = DMSetOptionsPrefix(user.dmpgrid,"pgrid_");CHKERRQ(ierr);
ierr = DMCompositeAddDM(user.dmpgrid,user.dmgen);CHKERRQ(ierr);
ierr = DMCompositeAddDM(user.dmpgrid,user.dmnet);CHKERRQ(ierr);
/* Create TAO solver and set desired solution method */
ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr);
ierr = TaoSetType(tao,TAOBLMVM);CHKERRQ(ierr);
/*
Optimization starts
*/
/* Set initial solution guess */
ierr = VecCreateSeq(PETSC_COMM_WORLD,3,&p);CHKERRQ(ierr);
ierr = VecGetArray(p,&x_ptr);CHKERRQ(ierr);
x_ptr[0] = PG[0]; x_ptr[1] = PG[1]; x_ptr[2] = PG[2];
ierr = VecRestoreArray(p,&x_ptr);CHKERRQ(ierr);
ierr = TaoSetInitialVector(tao,p);CHKERRQ(ierr);
/* Set routine for function and gradient evaluation */
ierr = TaoSetObjectiveRoutine(tao,FormFunction,(void *)&user);CHKERRQ(ierr);
ierr = TaoSetGradientRoutine(tao,TaoDefaultComputeGradient,(void *)&user);CHKERRQ(ierr);
/* Set bounds for the optimization */
ierr = VecDuplicate(p,&lowerb);CHKERRQ(ierr);
ierr = VecDuplicate(p,&upperb);CHKERRQ(ierr);
ierr = VecGetArray(lowerb,&x_ptr);CHKERRQ(ierr);
x_ptr[0] = 0.5; x_ptr[1] = 0.5; x_ptr[2] = 0.5;
ierr = VecRestoreArray(lowerb,&x_ptr);CHKERRQ(ierr);
ierr = VecGetArray(upperb,&x_ptr);CHKERRQ(ierr);
x_ptr[0] = 2.0; x_ptr[1] = 2.0; x_ptr[2] = 2.0;
ierr = VecRestoreArray(upperb,&x_ptr);CHKERRQ(ierr);
ierr = TaoSetVariableBounds(tao,lowerb,upperb);
/* Check for any TAO command line options */
ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);
ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr);
if (ksp) {
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
}
/* SOLVE THE APPLICATION */
ierr = TaoSolve(tao); CHKERRQ(ierr);
/* Get information on termination */
ierr = TaoGetConvergedReason(tao,&reason);CHKERRQ(ierr);
if (reason <= 0){
ierr=PetscPrintf(MPI_COMM_WORLD, "Try another method! \n");CHKERRQ(ierr);
}
ierr = VecView(p,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/* Free TAO data structures */
ierr = TaoDestroy(&tao);CHKERRQ(ierr);
ierr = VecDestroy(&user.vec_q);CHKERRQ(ierr);
ierr = VecDestroy(&lowerb);CHKERRQ(ierr);
ierr = VecDestroy(&upperb);CHKERRQ(ierr);
ierr = VecDestroy(&p);CHKERRQ(ierr);
ierr = DMDestroy(&user.dmgen);CHKERRQ(ierr);
ierr = DMDestroy(&user.dmnet);CHKERRQ(ierr);
ierr = DMDestroy(&user.dmpgrid);CHKERRQ(ierr);
ierr = ISDestroy(&user.is_diff);CHKERRQ(ierr);
ierr = ISDestroy(&user.is_alg);CHKERRQ(ierr);
ierr = PetscFinalize();
return(0);
}示例7: PCMGSetType
/*@C
PCMGSetLevels - Sets the number of levels to use with MG.
Must be called before any other MG routine.
Logically Collective on PC
Input Parameters:
+ pc - the preconditioner context
. levels - the number of levels
- comms - optional communicators for each level; this is to allow solving the coarser problems
on smaller sets of processors. Use NULL_OBJECT for default in Fortran
Level: intermediate
Notes:
If the number of levels is one then the multigrid uses the -mg_levels prefix
for setting the level options rather than the -mg_coarse prefix.
.keywords: MG, set, levels, multigrid
.seealso: PCMGSetType(), PCMGGetLevels()
@*/
PetscErrorCode PCMGSetLevels(PC pc,PetscInt levels,MPI_Comm *comms)
{
PetscErrorCode ierr;
PC_MG *mg = (PC_MG*)pc->data;
MPI_Comm comm;
PC_MG_Levels **mglevels = mg->levels;
PCMGType mgtype = mg->am;
PetscInt mgctype = (PetscInt) PC_MG_CYCLE_V;
PetscInt i;
PetscMPIInt size;
const char *prefix;
PC ipc;
PetscInt n;
PetscFunctionBegin;
PetscValidHeaderSpecific(pc,PC_CLASSID,1);
PetscValidLogicalCollectiveInt(pc,levels,2);
ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
if (mg->nlevels == levels) PetscFunctionReturn(0);
if (mglevels) {
mgctype = mglevels[0]->cycles;
/* changing the number of levels so free up the previous stuff */
ierr = PCReset_MG(pc);CHKERRQ(ierr);
n = mglevels[0]->levels;
for (i=0; i<n; i++) {
if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
ierr = KSPDestroy(&mglevels[i]->smoothd);CHKERRQ(ierr);
}
ierr = KSPDestroy(&mglevels[i]->smoothu);CHKERRQ(ierr);
ierr = PetscFree(mglevels[i]);CHKERRQ(ierr);
}
ierr = PetscFree(mg->levels);CHKERRQ(ierr);
}
mg->nlevels = levels;
ierr = PetscMalloc1(levels,&mglevels);CHKERRQ(ierr);
ierr = PetscLogObjectMemory((PetscObject)pc,levels*(sizeof(PC_MG*)));CHKERRQ(ierr);
ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
mg->stageApply = 0;
for (i=0; i<levels; i++) {
ierr = PetscNewLog(pc,&mglevels[i]);CHKERRQ(ierr);
mglevels[i]->level = i;
mglevels[i]->levels = levels;
mglevels[i]->cycles = mgctype;
mg->default_smoothu = 2;
mg->default_smoothd = 2;
mglevels[i]->eventsmoothsetup = 0;
mglevels[i]->eventsmoothsolve = 0;
mglevels[i]->eventresidual = 0;
mglevels[i]->eventinterprestrict = 0;
if (comms) comm = comms[i];
ierr = KSPCreate(comm,&mglevels[i]->smoothd);CHKERRQ(ierr);
ierr = KSPSetErrorIfNotConverged(mglevels[i]->smoothd,pc->erroriffailure);CHKERRQ(ierr);
ierr = PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd,(PetscObject)pc,levels-i);CHKERRQ(ierr);
ierr = KSPSetOptionsPrefix(mglevels[i]->smoothd,prefix);CHKERRQ(ierr);
ierr = PetscObjectComposedDataSetInt((PetscObject) mglevels[i]->smoothd, PetscMGLevelId, mglevels[i]->level);CHKERRQ(ierr);
if (i || levels == 1) {
char tprefix[128];
ierr = KSPSetType(mglevels[i]->smoothd,KSPCHEBYSHEV);CHKERRQ(ierr);
ierr = KSPSetConvergenceTest(mglevels[i]->smoothd,KSPConvergedSkip,NULL,NULL);CHKERRQ(ierr);
ierr = KSPSetNormType(mglevels[i]->smoothd,KSP_NORM_NONE);CHKERRQ(ierr);
ierr = KSPGetPC(mglevels[i]->smoothd,&ipc);CHKERRQ(ierr);
ierr = PCSetType(ipc,PCSOR);CHKERRQ(ierr);
ierr = KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, mg->default_smoothd);CHKERRQ(ierr);
sprintf(tprefix,"mg_levels_%d_",(int)i);
ierr = KSPAppendOptionsPrefix(mglevels[i]->smoothd,tprefix);CHKERRQ(ierr);
} else {
ierr = KSPAppendOptionsPrefix(mglevels[0]->smoothd,"mg_coarse_");CHKERRQ(ierr);
/* coarse solve is (redundant) LU by default; set shifttype NONZERO to avoid annoying zero-pivot in LU preconditioner */
ierr = KSPSetType(mglevels[0]->smoothd,KSPPREONLY);CHKERRQ(ierr);
//.........这里部分代码省略.........
示例8: output
//.........这里部分代码省略.........
ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr);
#if defined(PETSC_HAVE_SAWS)
ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);CHKERRQ(ierr);
#endif
if (iascii) {
ierr = PetscObjectPrintClassNamePrefixType((PetscObject)ksp,viewer);CHKERRQ(ierr);
if (ksp->ops->view) {
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
ierr = (*ksp->ops->view)(ksp,viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
}
if (ksp->guess_zero) {
ierr = PetscViewerASCIIPrintf(viewer," maximum iterations=%D, initial guess is zero\n",ksp->max_it);CHKERRQ(ierr);
} else {
ierr = PetscViewerASCIIPrintf(viewer," maximum iterations=%D\n", ksp->max_it);CHKERRQ(ierr);
}
if (ksp->guess_knoll) {ierr = PetscViewerASCIIPrintf(viewer," using preconditioner applied to right hand side for initial guess\n");CHKERRQ(ierr);}
ierr = PetscViewerASCIIPrintf(viewer," tolerances: relative=%g, absolute=%g, divergence=%g\n",(double)ksp->rtol,(double)ksp->abstol,(double)ksp->divtol);CHKERRQ(ierr);
if (ksp->pc_side == PC_RIGHT) {
ierr = PetscViewerASCIIPrintf(viewer," right preconditioning\n");CHKERRQ(ierr);
} else if (ksp->pc_side == PC_SYMMETRIC) {
ierr = PetscViewerASCIIPrintf(viewer," symmetric preconditioning\n");CHKERRQ(ierr);
} else {
ierr = PetscViewerASCIIPrintf(viewer," left preconditioning\n");CHKERRQ(ierr);
}
if (ksp->guess) {ierr = PetscViewerASCIIPrintf(viewer," using Fischers initial guess method %D with size %D\n",ksp->guess->method,ksp->guess->maxl);CHKERRQ(ierr);}
if (ksp->dscale) {ierr = PetscViewerASCIIPrintf(viewer," diagonally scaled system\n");CHKERRQ(ierr);}
if (!ksp->guess_zero) {ierr = PetscViewerASCIIPrintf(viewer," using nonzero initial guess\n");CHKERRQ(ierr);}
ierr = PetscViewerASCIIPrintf(viewer," using %s norm type for convergence test\n",KSPNormTypes[ksp->normtype]);CHKERRQ(ierr);
} else if (isbinary) {
PetscInt classid = KSP_FILE_CLASSID;
MPI_Comm comm;
PetscMPIInt rank;
char type[256];
ierr = PetscObjectGetComm((PetscObject)ksp,&comm);CHKERRQ(ierr);
ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
if (!rank) {
ierr = PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr);
ierr = PetscStrncpy(type,((PetscObject)ksp)->type_name,256);CHKERRQ(ierr);
ierr = PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);CHKERRQ(ierr);
}
if (ksp->ops->view) {
ierr = (*ksp->ops->view)(ksp,viewer);CHKERRQ(ierr);
}
} else if (isdraw) {
PetscDraw draw;
char str[36];
PetscReal x,y,bottom,h;
PetscBool flg;
ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr);
ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);CHKERRQ(ierr);
if (!flg) {
ierr = PetscStrcpy(str,"KSP: ");CHKERRQ(ierr);
ierr = PetscStrcat(str,((PetscObject)ksp)->type_name);CHKERRQ(ierr);
ierr = PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr);
bottom = y - h;
} else {
bottom = y;
}
ierr = PetscDrawPushCurrentPoint(draw,x,bottom);CHKERRQ(ierr);
#if defined(PETSC_HAVE_SAWS)
} else if (issaws) {
PetscMPIInt rank;
const char *name;
ierr = PetscObjectGetName((PetscObject)ksp,&name);CHKERRQ(ierr);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
if (!((PetscObject)ksp)->amsmem && !rank) {
char dir[1024];
ierr = PetscObjectViewSAWs((PetscObject)ksp,viewer);CHKERRQ(ierr);
ierr = PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/its",name);CHKERRQ(ierr);
PetscStackCallSAWs(SAWs_Register,(dir,&ksp->its,1,SAWs_READ,SAWs_INT));
if (!ksp->res_hist) {
ierr = KSPSetResidualHistory(ksp,NULL,PETSC_DECIDE,PETSC_TRUE);CHKERRQ(ierr);
}
ierr = PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/res_hist",name);CHKERRQ(ierr);
PetscStackCallSAWs(SAWs_Register,(dir,ksp->res_hist,10,SAWs_READ,SAWs_DOUBLE));
}
#endif
} else if (ksp->ops->view) {
ierr = (*ksp->ops->view)(ksp,viewer);CHKERRQ(ierr);
}
if (!ksp->skippcsetfromoptions) {
if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);}
ierr = PCView(ksp->pc,viewer);CHKERRQ(ierr);
}
if (isdraw) {
PetscDraw draw;
ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}示例9: port_lsd_bfbt
PetscErrorCode port_lsd_bfbt(void)
{
Mat A;
Vec x,b;
KSP ksp_A;
PC pc_A;
IS isu,isp;
PetscErrorCode ierr;
PetscFunctionBeginUser;
ierr = LoadTestMatrices(&A,&x,&b,&isu,&isp);CHKERRQ(ierr);
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp_A);CHKERRQ(ierr);
ierr = KSPSetOptionsPrefix(ksp_A,"fc_");CHKERRQ(ierr);
ierr = KSPSetOperators(ksp_A,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = KSPGetPC(ksp_A,&pc_A);CHKERRQ(ierr);
ierr = PCSetType(pc_A,PCFIELDSPLIT);CHKERRQ(ierr);
ierr = PCFieldSplitSetBlockSize(pc_A,2);CHKERRQ(ierr);
ierr = PCFieldSplitSetIS(pc_A,"velocity",isu);CHKERRQ(ierr);
ierr = PCFieldSplitSetIS(pc_A,"pressure",isp);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp_A);CHKERRQ(ierr);
ierr = KSPSolve(ksp_A,b,x);CHKERRQ(ierr);
/* Pull u,p out of x */
{
PetscInt loc;
PetscReal max,norm;
PetscScalar sum;
Vec uvec,pvec;
VecScatter uscat,pscat;
Mat A11,A22;
/* grab matrices and create the compatable u,p vectors */
ierr = MatGetSubMatrix(A,isu,isu,MAT_INITIAL_MATRIX,&A11);CHKERRQ(ierr);
ierr = MatGetSubMatrix(A,isp,isp,MAT_INITIAL_MATRIX,&A22);CHKERRQ(ierr);
ierr = MatGetVecs(A11,&uvec,PETSC_NULL);CHKERRQ(ierr);
ierr = MatGetVecs(A22,&pvec,PETSC_NULL);CHKERRQ(ierr);
/* perform the scatter from x -> (u,p) */
ierr = VecScatterCreate(x,isu,uvec,PETSC_NULL,&uscat);CHKERRQ(ierr);
ierr = VecScatterBegin(uscat,x,uvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecScatterEnd(uscat,x,uvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecScatterCreate(x,isp,pvec,PETSC_NULL,&pscat);CHKERRQ(ierr);
ierr = VecScatterBegin(pscat,x,pvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecScatterEnd(pscat,x,pvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD,"-- vector vector values --\n");
ierr = VecMin(uvec,&loc,&max);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Min(u) = %1.6F [loc=%d]\n",max,loc);CHKERRQ(ierr);
ierr = VecMax(uvec,&loc,&max);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Max(u) = %1.6F [loc=%d]\n",max,loc);CHKERRQ(ierr);
ierr = VecNorm(uvec,NORM_2,&norm);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Norm(u) = %1.6F \n",norm);CHKERRQ(ierr);
ierr = VecSum(uvec,&sum);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Sum(u) = %1.6F \n",PetscRealPart(sum));CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD,"-- pressure vector values --\n");
ierr = VecMin(pvec,&loc,&max);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Min(p) = %1.6F [loc=%d]\n",max,loc);CHKERRQ(ierr);
ierr = VecMax(pvec,&loc,&max);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Max(p) = %1.6F [loc=%d]\n",max,loc);CHKERRQ(ierr);
ierr = VecNorm(pvec,NORM_2,&norm);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Norm(p) = %1.6F \n",norm);CHKERRQ(ierr);
ierr = VecSum(pvec,&sum);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Sum(p) = %1.6F \n",PetscRealPart(sum));CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD,"-- Full vector values --\n");
ierr = VecMin(x,&loc,&max);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Min(u,p) = %1.6F [loc=%d]\n",max,loc);CHKERRQ(ierr);
ierr = VecMax(x,&loc,&max);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Max(u,p) = %1.6F [loc=%d]\n",max,loc);CHKERRQ(ierr);
ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Norm(u,p) = %1.6F \n",norm);CHKERRQ(ierr);
ierr = VecSum(x,&sum);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Sum(u,p) = %1.6F \n",PetscRealPart(sum));CHKERRQ(ierr);
ierr = VecScatterDestroy(&uscat);CHKERRQ(ierr);
ierr = VecScatterDestroy(&pscat);CHKERRQ(ierr);
ierr = VecDestroy(&uvec);CHKERRQ(ierr);
ierr = VecDestroy(&pvec);CHKERRQ(ierr);
ierr = MatDestroy(&A11);CHKERRQ(ierr);
ierr = MatDestroy(&A22);CHKERRQ(ierr);
}
ierr = KSPDestroy(&ksp_A);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr);
ierr = ISDestroy(&isu);CHKERRQ(ierr);
ierr = ISDestroy(&isp);CHKERRQ(ierr);
PetscFunctionReturn(0);
}示例10: KSPSolve_FBCGSR
static PetscErrorCode KSPSolve_FBCGSR(KSP ksp)
{
PetscErrorCode ierr;
PetscInt i,j,N;
PetscScalar tau,sigma,alpha,omega,beta;
PetscReal rho;
PetscScalar xi1,xi2,xi3,xi4;
Vec X,B,P,P2,RP,R,V,S,T,S2;
PetscScalar *PETSC_RESTRICT rp, *PETSC_RESTRICT r, *PETSC_RESTRICT p;
PetscScalar *PETSC_RESTRICT v, *PETSC_RESTRICT s, *PETSC_RESTRICT t, *PETSC_RESTRICT s2;
PetscScalar insums[4],outsums[4];
KSP_BCGS *bcgs = (KSP_BCGS*)ksp->data;
PC pc;
Mat mat;
PetscFunctionBegin;
if (!ksp->vec_rhs->petscnative) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Only coded for PETSc vectors");
ierr = VecGetLocalSize(ksp->vec_sol,&N);CHKERRQ(ierr);
X = ksp->vec_sol;
B = ksp->vec_rhs;
P2 = ksp->work[0];
/* The followings are involved in modified inner product calculations and vector updates */
RP = ksp->work[1]; ierr = VecGetArray(RP,(PetscScalar**)&rp);CHKERRQ(ierr); ierr = VecRestoreArray(RP,NULL);CHKERRQ(ierr);
R = ksp->work[2]; ierr = VecGetArray(R,(PetscScalar**)&r);CHKERRQ(ierr); ierr = VecRestoreArray(R,NULL);CHKERRQ(ierr);
P = ksp->work[3]; ierr = VecGetArray(P,(PetscScalar**)&p);CHKERRQ(ierr); ierr = VecRestoreArray(P,NULL);CHKERRQ(ierr);
V = ksp->work[4]; ierr = VecGetArray(V,(PetscScalar**)&v);CHKERRQ(ierr); ierr = VecRestoreArray(V,NULL);CHKERRQ(ierr);
S = ksp->work[5]; ierr = VecGetArray(S,(PetscScalar**)&s);CHKERRQ(ierr); ierr = VecRestoreArray(S,NULL);CHKERRQ(ierr);
T = ksp->work[6]; ierr = VecGetArray(T,(PetscScalar**)&t);CHKERRQ(ierr); ierr = VecRestoreArray(T,NULL);CHKERRQ(ierr);
S2 = ksp->work[7]; ierr = VecGetArray(S2,(PetscScalar**)&s2);CHKERRQ(ierr); ierr = VecRestoreArray(S2,NULL);CHKERRQ(ierr);
/* Only supports right preconditioning */
if (ksp->pc_side != PC_RIGHT) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSP fbcgsr does not support %s",PCSides[ksp->pc_side]);
if (!ksp->guess_zero) {
if (!bcgs->guess) {
ierr = VecDuplicate(X,&bcgs->guess);CHKERRQ(ierr);
}
ierr = VecCopy(X,bcgs->guess);CHKERRQ(ierr);
} else {
ierr = VecSet(X,0.0);CHKERRQ(ierr);
}
/* Compute initial residual */
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetUp(pc);CHKERRQ(ierr);
ierr = PCGetOperators(pc,&mat,NULL);CHKERRQ(ierr);
if (!ksp->guess_zero) {
ierr = KSP_MatMult(ksp,mat,X,P2);CHKERRQ(ierr); /* P2 is used as temporary storage */
ierr = VecCopy(B,R);CHKERRQ(ierr);
ierr = VecAXPY(R,-1.0,P2);CHKERRQ(ierr);
} else {
ierr = VecCopy(B,R);CHKERRQ(ierr);
}
/* Test for nothing to do */
ierr = VecNorm(R,NORM_2,&rho);CHKERRQ(ierr);
ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
ksp->its = 0;
ksp->rnorm = rho;
ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
ierr = KSPLogResidualHistory(ksp,rho);CHKERRQ(ierr);
ierr = KSPMonitor(ksp,0,rho);CHKERRQ(ierr);
ierr = (*ksp->converged)(ksp,0,rho,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
if (ksp->reason) PetscFunctionReturn(0);
/* Initialize iterates */
ierr = VecCopy(R,RP);CHKERRQ(ierr); /* rp <- r */
ierr = VecCopy(R,P);CHKERRQ(ierr); /* p <- r */
/* Big loop */
for (i=0; i<ksp->max_it; i++) {
/* matmult and pc */
ierr = KSP_PCApply(ksp,P,P2);CHKERRQ(ierr); /* p2 <- K p */
ierr = KSP_MatMult(ksp,mat,P2,V);CHKERRQ(ierr); /* v <- A p2 */
/* inner prodcuts */
if (i==0) {
tau = rho*rho;
ierr = VecDot(V,RP,&sigma);CHKERRQ(ierr); /* sigma <- (v,rp) */
} else {
ierr = PetscLogEventBegin(VEC_ReduceArithmetic,0,0,0,0);CHKERRQ(ierr);
tau = sigma = 0.0;
for (j=0; j<N; j++) {
tau += r[j]*rp[j]; /* tau <- (r,rp) */
sigma += v[j]*rp[j]; /* sigma <- (v,rp) */
}
ierr = PetscLogFlops(4.0*N);CHKERRQ(ierr);
ierr = PetscLogEventEnd(VEC_ReduceArithmetic,0,0,0,0);CHKERRQ(ierr);
insums[0] = tau;
insums[1] = sigma;
ierr = PetscLogEventBegin(VEC_ReduceCommunication,0,0,0,0);CHKERRQ(ierr);
ierr = MPIU_Allreduce(insums,outsums,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ksp));CHKERRQ(ierr);
ierr = PetscLogEventEnd(VEC_ReduceCommunication,0,0,0,0);CHKERRQ(ierr);
tau = outsums[0];
sigma = outsums[1];
}
/* scalar update */
//.........这里部分代码省略.........
示例11: main
int main(int argc,char **argv)
{
SNES snes; /* nonlinear solver context */
KSP ksp; /* linear solver context */
PC pc; /* preconditioner context */
Vec x,r; /* solution, residual vectors */
Mat J; /* Jacobian matrix */
PetscErrorCode ierr;
PetscInt its;
PetscMPIInt size;
PetscScalar pfive = .5,*xx;
PetscBool flg;
PetscInitialize(&argc,&argv,(char*)0,help);
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create nonlinear solver context
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create matrix and vector data structures; set corresponding routines
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Create vectors for solution and nonlinear function
*/
ierr = VecCreateSeq(PETSC_COMM_SELF,2,&x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&r);CHKERRQ(ierr);
/*
Create Jacobian matrix data structure
*/
ierr = MatCreate(PETSC_COMM_SELF,&J);CHKERRQ(ierr);
ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr);
ierr = MatSetFromOptions(J);CHKERRQ(ierr);
ierr = PetscOptionsHasName(NULL,"-hard",&flg);CHKERRQ(ierr);
if (!flg) {
/*
Set function evaluation routine and vector.
*/
ierr = SNESSetFunction(snes,r,FormFunction1,NULL);CHKERRQ(ierr);
/*
Set Jacobian matrix data structure and Jacobian evaluation routine
*/
ierr = SNESSetJacobian(snes,J,J,FormJacobian1,NULL);CHKERRQ(ierr);
} else {
ierr = SNESSetFunction(snes,r,FormFunction2,NULL);CHKERRQ(ierr);
ierr = SNESSetJacobian(snes,J,J,FormJacobian2,NULL);CHKERRQ(ierr);
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Customize nonlinear solver; set runtime options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Set linear solver defaults for this problem. By extracting the
KSP, KSP, and PC contexts from the SNES context, we can then
directly call any KSP, KSP, and PC routines to set various options.
*/
ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
ierr = KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);CHKERRQ(ierr);
/*
Set SNES/KSP/KSP/PC runtime options, e.g.,
-snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
These options will override those specified above as long as
SNESSetFromOptions() is called _after_ any other customization
routines.
*/
ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Evaluate initial guess; then solve nonlinear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
if (!flg) {
ierr = VecSet(x,pfive);CHKERRQ(ierr);
} else {
ierr = VecGetArray(x,&xx);CHKERRQ(ierr);
xx[0] = 2.0; xx[1] = 3.0;
ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr);
}
/*
Note: The user should initialize the vector, x, with the initial guess
for the nonlinear solver prior to calling SNESSolve(). In particular,
to employ an initial guess of zero, the user should explicitly set
this vector to zero by calling VecSet().
*/
ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr);
ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
if (flg) {
Vec f;
//.........这里部分代码省略.........
示例12: _solve_routine
/**
* Routine for solving a linear equation with PETSc.
*
* Note that this routine calls MPI teardown routines, so it should
* not be called in the main process unless you want to break MPI for
* the remainder of the process lifetime.
*
*/
PetscErrorCode _solve_routine(JNIEnv *env, jint n, jintArray *index,
jobjectArray *diagonals, jobjectArray *options,
double *solution, jdoubleArray *rhs) {
PetscErrorCode ierr;
KSP ksp;
PC pc;
Mat A;
Vec b;
Vec x;
PetscInitialize(0, NULL, (char*) NULL, NULL);
// Set up matrix A in Ax = b
int num_diags = (*env) -> GetArrayLength(env, *index);
ierr = MatCreateSeqAIJ(PETSC_COMM_WORLD, n, n, num_diags, NULL,
&A); CHKERRQ(ierr);
ierr = MatSetUp(A); CHKERRQ(ierr);
ierr = _fill_matrix(env, &A, index, diagonals);
CHKERRQ(ierr);
ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
// Set up vectors b and x in Ax = b
ierr = VecCreate(PETSC_COMM_WORLD, &b); CHKERRQ(ierr);
ierr = VecSetSizes(b, PETSC_DECIDE, n); CHKERRQ(ierr);
ierr = VecSetType(b, VECSEQ);
ierr = VecDuplicate(b, &x);
ierr = _fill_vector(env, &b, rhs); CHKERRQ(ierr);
ierr = VecAssemblyBegin(b); CHKERRQ(ierr);
ierr = VecAssemblyEnd(b); CHKERRQ(ierr);
// Set up solver
ierr = KSPCreate(PETSC_COMM_WORLD, &ksp); CHKERRQ(ierr);
ierr = KSPSetOperators(ksp, A, A); CHKERRQ(ierr);
// Solver options
ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr);
for(int i = 0; i < (*env)->GetArrayLength(env, *options); i++) {
jobject option = (*env)->GetObjectArrayElement(env, *options, i);
jstring key = (*env)->GetObjectArrayElement(env, option, 0);
jstring value = (*env)->GetObjectArrayElement(env, option, 1);
const jchar *keyChars = (*env)->GetStringChars(env, key, NULL);
const jchar *valChars = (*env)->GetStringChars(env, value, NULL);
PetscOptionsSetValue((char*) keyChars, (char*) valChars);
(*env)->ReleaseStringChars(env, key, keyChars);
(*env)->ReleaseStringChars(env, value, valChars);
}
ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
ierr = KSPSetUp(ksp); CHKERRQ(ierr);
// Solve
ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
// Copy solution into JNI allocated memory
ierr = _pvec_to_array(&x, solution); CHKERRQ(ierr);
// Clean up and return
ierr = MatDestroy(&A); CHKERRQ(ierr);
ierr = VecDestroy(&b); CHKERRQ(ierr);
ierr = VecDestroy(&x); CHKERRQ(ierr);
ierr = KSPDestroy(&ksp);
PetscFinalize();
return ierr;
}示例13: main
int main(int argc, char *argv[])
{
// Initialize libMesh
libMesh::LibMeshInit init(argc, argv);
libMesh::Parallel::Communicator& WorldComm = init.comm();
libMesh::PetscMatrix<libMesh::Number> matrix_A(WorldComm);
matrix_A.init(4,4,4,4);
matrix_A.set(0,0,1.);
// matrix_A.set(0,1,2.);
// matrix_A.set(0,2,3.);
// matrix_A.set(0,3,4.);
// matrix_A.set(1,0,2.);
matrix_A.set(1,1,5.);
// matrix_A.set(1,2,3.);
// matrix_A.set(1,3,7.);
// matrix_A.set(2,0,3.);
// matrix_A.set(2,1,3.);
matrix_A.set(2,2,9.);
// matrix_A.set(2,3,6.);
// matrix_A.set(3,0,4.);
// matrix_A.set(3,1,7.);
// matrix_A.set(3,2,6.);
matrix_A.set(3,3,1.);
matrix_A.close();
Mat dummy_inv_A;
MatCreate(PETSC_COMM_WORLD,&dummy_inv_A);
MatSetType(dummy_inv_A,MATMPIAIJ);
MatSetSizes(dummy_inv_A,PETSC_DECIDE,PETSC_DECIDE,4,4);
MatMPIAIJSetPreallocation(dummy_inv_A,2,NULL,0,NULL);
MatSetUp(dummy_inv_A);
// Dummy matrices
// Mat dummy_A, dummy_inv_A;
//
//
libMesh::PetscVector<libMesh::Number> vector_unity(WorldComm,4,4);
libMesh::PetscVector<libMesh::Number> vector_dummy_answer(WorldComm,4,4);
VecSet(vector_unity.vec(),1);
vector_unity.close();
VecSet(vector_dummy_answer.vec(),0);
vector_dummy_answer.close();
// Solver
// libMesh::PetscLinearSolver<libMesh::Number> KSP_dummy_solver(WorldComm);
// KSP_dummy_solver.init(&matrix_A);
// KSPSetOperators(KSP_dummy_solver.ksp(),matrix_A.mat(),NULL);
KSP ksp;
PC pc;
KSPCreate(PETSC_COMM_WORLD,&ksp);
KSPSetOperators(ksp, matrix_A.mat(), matrix_A.mat());
KSPGetPC(ksp,&pc);
PCSetFromOptions(pc);
PCType dummy_type;
PCGetType(pc,&dummy_type);
std::cout << std::endl << dummy_type << std::endl << std::endl;
// PCSetType(pc,PCSPAI);
// PCHYPRESetType(pc,"parasails");
KSPSetUp(ksp);
KSPSolve(ksp,vector_unity.vec(),vector_dummy_answer.vec());
PCComputeExplicitOperator(pc,&dummy_inv_A);
// KSPGetOperators(KSP_dummy_solver.ksp(),&dummy_A,&dummy_inv_A);
libMesh::PetscMatrix<libMesh::Number> matrix_invA(dummy_inv_A,WorldComm);
matrix_invA.close();
//
// // KSP_dummy_solver.solve(matrix_A,vector_dummy_answer,vector_unity,1E-5,10000);
//
// vector_dummy_answer.print_matlab();
//
// libMesh::PetscMatrix<libMesh::Number> product_mat(WorldComm);
matrix_A.print_matlab();
matrix_invA.print_matlab();
vector_dummy_answer.print_matlab();
return 0;
}示例14: main
int main(int argc, char **argv){
PetscErrorCode ierr;
int nx = 63, ny = 63;
DM dm;
PetscBool flg;
Mat A;
Vec u, b;
KSP solver;
PC pc;
double norm;
PetscInt stage;
ierr = PetscInitialize(&argc, &argv, NULL, NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL, "-nx", &nx, PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL, "-ny", &ny, PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsHasName(PETSC_NULL, "-assemble", &flg);CHKERRQ(ierr);
ierr = PetscLogStageRegister("preparing",&stage); CHKERRQ(ierr);
ierr = PetscLogStagePush(stage);CHKERRQ(ierr);
ierr = PetscLogStageRegister("Domain creation",&stage); CHKERRQ(ierr);
ierr = PetscLogStagePush(stage);CHKERRQ(ierr);
ierr = createDomain(&dm, nx, ny);CHKERRQ(ierr);
ierr = PetscLogStagePop();CHKERRQ(ierr);
ierr = PetscLogStageRegister("matrix creation",&stage); CHKERRQ(ierr);
ierr = PetscLogStagePush(stage);CHKERRQ(ierr);
ierr = createMat(dm, &A, flg);CHKERRQ(ierr);
ierr = PetscLogStagePop();CHKERRQ(ierr);
ierr = PetscLogStageRegister("Vector creation",&stage); CHKERRQ(ierr);
ierr = PetscLogStagePush(stage);CHKERRQ(ierr);
ierr = PetscLogStagePop();CHKERRQ(ierr);
ierr = DMCreateGlobalVector(dm, &b);CHKERRQ(ierr);
ierr = VecDuplicate(b, &u);CHKERRQ(ierr);
ierr = PetscLogStageRegister("Domain initialisation",&stage); CHKERRQ(ierr);
ierr = PetscLogStagePush(stage);CHKERRQ(ierr);
ierr = init2d(dm, b);CHKERRQ(ierr);
ierr = PetscLogStagePop();CHKERRQ(ierr);
ierr = PetscLogStageRegister("solver creation",&stage); CHKERRQ(ierr);
ierr = PetscLogStagePush(stage);CHKERRQ(ierr);
ierr = KSPCreate(PETSC_COMM_WORLD, &solver);CHKERRQ(ierr);
ierr = KSPSetOptionsPrefix(solver, "poisson_");CHKERRQ(ierr);
ierr = KSPSetOperators(solver, A, A, DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = KSPSetType(solver, KSPCG);
ierr = KSPGetPC(solver, &pc);CHKERRQ(ierr);
ierr = PCSetType(pc, PCNONE);CHKERRQ(ierr);
ierr = KSPSetFromOptions(solver);CHKERRQ(ierr);
ierr = PetscLogStagePop();CHKERRQ(ierr);
ierr = PetscLogStagePop();CHKERRQ(ierr);
ierr = PetscLogStageRegister("Solving",&stage); CHKERRQ(ierr);
ierr = PetscLogStagePush(stage);CHKERRQ(ierr);
ierr = KSPSolve(solver, b, u);CHKERRQ(ierr);
ierr = PetscLogStagePop();CHKERRQ(ierr);
//ierr = VecView(u, PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
//sleep(10);
VecDestroy(&u);
VecDestroy(&b);
MatDestroy(&A);
DMDestroy(&dm);
KSPDestroy(&solver);
ierr = PetscFinalize();
return 0;
}示例15: main
//.........这里部分代码省略.........
ierr = PetscOptionsHasName(NULL,"-ts_fd",&flg);CHKERRQ(ierr);
if (!flg) {
ierr = TSSetRHSJacobian(ts,J,J,RHSJacobian,&data);CHKERRQ(ierr);
} else {
ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
ierr = PetscOptionsHasName(NULL,"-fd_color",&fd_jacobian_coloring);CHKERRQ(ierr);
if (fd_jacobian_coloring) { /* Use finite differences with coloring */
/* Get data structure of J */
PetscBool pc_diagonal;
ierr = PetscOptionsHasName(NULL,"-pc_diagonal",&pc_diagonal);CHKERRQ(ierr);
if (pc_diagonal) { /* the preconditioner of J is a diagonal matrix */
PetscInt rstart,rend,i;
PetscScalar zero=0.0;
ierr = MatGetOwnershipRange(J,&rstart,&rend);CHKERRQ(ierr);
for (i=rstart; i<rend; i++) {
ierr = MatSetValues(J,1,&i,1,&i,&zero,INSERT_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
} else {
/* Fill the structure using the expensive SNESComputeJacobianDefault. Temporarily set up the TS so we can call this function */
ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr);
ierr = TSSetUp(ts);CHKERRQ(ierr);
ierr = SNESComputeJacobianDefault(snes,x,J,J,ts);CHKERRQ(ierr);
}
/* create coloring context */
ierr = MatColoringCreate(J,&mc);CHKERRQ(ierr);
ierr = MatColoringSetType(mc,MATCOLORINGSL);CHKERRQ(ierr);
ierr = MatColoringSetFromOptions(mc);CHKERRQ(ierr);
ierr = MatColoringApply(mc,&iscoloring);CHKERRQ(ierr);
ierr = MatColoringDestroy(&mc);CHKERRQ(ierr);
ierr = MatFDColoringCreate(J,iscoloring,&matfdcoloring);CHKERRQ(ierr);
ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts);CHKERRQ(ierr);
ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr);
ierr = MatFDColoringSetUp(J,iscoloring,matfdcoloring);CHKERRQ(ierr);
ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);CHKERRQ(ierr);
ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr);
} else { /* Use finite differences (slow) */
ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,NULL);CHKERRQ(ierr);
}
}
/* Pick up a Petsc preconditioner */
/* one can always set method or preconditioner during the run time */
ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
ierr = TSSetUp(ts);CHKERRQ(ierr);
/* Test TSSetPostStep() */
ierr = PetscOptionsHasName(NULL,"-test_PostStep",&flg);CHKERRQ(ierr);
if (flg) {
ierr = TSSetPostStep(ts,PostStep);CHKERRQ(ierr);
}
ierr = PetscOptionsGetInt(NULL,"-NOUT",&NOUT,NULL);CHKERRQ(ierr);
for (iout=1; iout<=NOUT; iout++) {
ierr = TSSetDuration(ts,time_steps,iout*ftime_original/NOUT);CHKERRQ(ierr);
ierr = TSSolve(ts,global);CHKERRQ(ierr);
ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);
ierr = TSSetInitialTimeStep(ts,ftime,dt);CHKERRQ(ierr);
}
/* Interpolate solution at tfinal */
ierr = TSGetSolution(ts,&global);CHKERRQ(ierr);
ierr = TSInterpolate(ts,ftime_original,global);CHKERRQ(ierr);
ierr = PetscOptionsHasName(NULL,"-matlab_view",&flg);CHKERRQ(ierr);
if (flg) { /* print solution into a MATLAB file */
ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile);CHKERRQ(ierr);
ierr = PetscViewerSetFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);
ierr = VecView(global,viewfile);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewfile);CHKERRQ(ierr);
}
/* display solver info for Sundials */
ierr = TSGetType(ts,&tstype);CHKERRQ(ierr);
ierr = PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&sundials);CHKERRQ(ierr);
if (sundials) {
ierr = PetscViewerStringOpen(PETSC_COMM_WORLD,tsinfo,120,&viewer);CHKERRQ(ierr);
ierr = TSView(ts,viewer);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
ierr = PetscViewerStringOpen(PETSC_COMM_WORLD,pcinfo,120,&viewer);CHKERRQ(ierr);
ierr = PCView(pc,viewer);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"%d Procs,%s TSType, %s Preconditioner\n",size,tsinfo,pcinfo);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
}
/* free the memories */
ierr = TSDestroy(&ts);CHKERRQ(ierr);
ierr = VecDestroy(&global);CHKERRQ(ierr);
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = MatDestroy(&J);CHKERRQ(ierr);
if (fd_jacobian_coloring) {ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr);}
ierr = PetscFinalize();
return 0;
}