本文整理汇总了C++中KSPSolve函数的典型用法代码示例。如果您正苦于以下问题:C++ KSPSolve函数的具体用法?C++ KSPSolve怎么用?C++ KSPSolve使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了KSPSolve函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: KSPSolve_SpecEst
static PetscErrorCode KSPSolve_SpecEst(KSP ksp)
{
PetscErrorCode ierr;
KSP_SpecEst *spec = (KSP_SpecEst*)ksp->data;
PetscFunctionBegin;
if (spec->current) {
ierr = KSPSolve(spec->kspcheap,ksp->vec_rhs,ksp->vec_sol);CHKERRQ(ierr);
ierr = KSPSpecEstPropagateUp(ksp,spec->kspcheap);CHKERRQ(ierr);
} else {
PetscInt i,its,neig;
PetscReal *real,*imag,rad = 0;
ierr = KSPSolve(spec->kspest,ksp->vec_rhs,ksp->vec_sol);CHKERRQ(ierr);
ierr = KSPSpecEstPropagateUp(ksp,spec->kspest);CHKERRQ(ierr);
ierr = KSPComputeExtremeSingularValues(spec->kspest,&spec->max,&spec->min);CHKERRQ(ierr);
ierr = KSPGetIterationNumber(spec->kspest,&its);CHKERRQ(ierr);
ierr = PetscMalloc2(its,PetscReal,&real,its,PetscReal,&imag);CHKERRQ(ierr);
ierr = KSPComputeEigenvalues(spec->kspest,its,real,imag,&neig);CHKERRQ(ierr);
for (i=0; i<neig; i++) {
/* We would really like to compute w (nominally 1/radius) to minimize |1-wB|. Empirically it
is better to compute rad = |1-B| than rad = |B|. There must be a cheap way to do better. */
rad = PetscMax(rad,PetscRealPart(PetscSqrtScalar((PetscScalar)(PetscSqr(real[i]-1.) + PetscSqr(imag[i])))));
}
ierr = PetscFree2(real,imag);CHKERRQ(ierr);
spec->radius = rad;
ierr = KSPChebyshevSetEigenvalues(spec->kspcheap,spec->max*spec->maxfactor,spec->min*spec->minfactor);CHKERRQ(ierr);
ierr = KSPRichardsonSetScale(spec->kspcheap,spec->richfactor/spec->radius);
ierr = PetscInfo3(ksp,"Estimated singular value min=%G max=%G, spectral radius=%G",spec->min,spec->max,spec->radius);CHKERRQ(ierr);
spec->current = PETSC_TRUE;
}
PetscFunctionReturn(0);
}
开发者ID:Kun-Qu,项目名称:petsc,代码行数:34,代码来源:specest.c
示例2: PCMGKCycle_Private
PetscErrorCode PCMGKCycle_Private(PC pc,PC_MG_Levels **mglevels)
{
PetscErrorCode ierr;
PetscInt i,l = mglevels[0]->levels;
PetscFunctionBegin;
/* restrict the RHS through all levels to coarsest. */
for (i=l-1; i>0; i--){
if (mglevels[i]->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
ierr = MatRestrict(mglevels[i]->restrct,mglevels[i]->b,mglevels[i-1]->b);CHKERRQ(ierr);
if (mglevels[i]->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
}
/* work our way up through the levels */
ierr = VecSet(mglevels[0]->x,0.0);CHKERRQ(ierr);
for (i=0; i<l-1; i++) {
if (mglevels[i]->eventsmoothsolve) {ierr = PetscLogEventBegin(mglevels[i]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
ierr = KSPSolve(mglevels[i]->smoothd,mglevels[i]->b,mglevels[i]->x);CHKERRQ(ierr);
if (mglevels[i]->eventsmoothsolve) {ierr = PetscLogEventEnd(mglevels[i]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
if (mglevels[i+1]->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels[i+1]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
ierr = MatInterpolate(mglevels[i+1]->interpolate,mglevels[i]->x,mglevels[i+1]->x);CHKERRQ(ierr);
if (mglevels[i+1]->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels[i+1]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
}
if (mglevels[l-1]->eventsmoothsolve) {ierr = PetscLogEventBegin(mglevels[l-1]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
ierr = KSPSolve(mglevels[l-1]->smoothd,mglevels[l-1]->b,mglevels[l-1]->x);CHKERRQ(ierr);
if (mglevels[l-1]->eventsmoothsolve) {ierr = PetscLogEventEnd(mglevels[l-1]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
PetscFunctionReturn(0);
}
开发者ID:Kun-Qu,项目名称:petsc,代码行数:29,代码来源:fmg.c
示例3: G
/*
S^{-1} = ( G^T G )^{-1} G^T K G ( G^T G )^{-1}
= A C A
S^{-T} = A^T (A C)^T
= A^T C^T A^T, but A = G^T G which is symmetric
= A C^T A
= A G^T ( G^T K )^T A
= A G^T K^T G A
*/
PetscErrorCode BSSCR_PCApplyTranspose_GtKG( PC pc, Vec x, Vec y )
{
PC_GtKG ctx = (PC_GtKG)pc->data;
KSP ksp;
Mat K, G;
Vec s,t,X;
PetscLogDouble t0,t1;
ksp = ctx->ksp;
K = ctx->K;
G = ctx->G;
s = ctx->s;
t = ctx->t;
X = ctx->X;
if (ctx->monitor_rhs_consistency) { BSSCRBSSCR_Lp_monitor_check_rhs_consistency(ksp,x,1); }
PetscGetTime(&t0);
KSPSolve( ksp, x, t ); /* t <- GtG_inv x */
PetscGetTime(&t1);
if (ctx->monitor_activated) { BSSCR_PCBFBTSubKSPMonitor(ksp,1,(t1-t0)); }
MatMult( G, t, s ); /* s <- G t */
MatMultTranspose( K, s, X ); /* X <- K^T s */
MatMultTranspose( G, X, t ); /* t <- Gt X */
if (ctx->monitor_rhs_consistency) { BSSCRBSSCR_Lp_monitor_check_rhs_consistency(ksp,t,2); }
PetscGetTime(&t0);
KSPSolve( ksp, t, y ); /* y <- GtG_inv t */
PetscGetTime(&t1);
if (ctx->monitor_activated) { BSSCR_PCBFBTSubKSPMonitor(ksp,2,(t1-t0)); }
PetscFunctionReturn(0);
}
开发者ID:AngelValverdePerez,项目名称:underworld2,代码行数:45,代码来源:pc_GtKG.c
示例4: BSSCR_BSSCR_PCApplyTranspose_GtKG_diagonal_scaling
PetscErrorCode BSSCR_BSSCR_PCApplyTranspose_GtKG_diagonal_scaling( PC pc, Vec x, Vec y )
{
PC_GtKG ctx = (PC_GtKG)pc->data;
KSP ksp;
Mat K, G;
Vec s,t,X,di;
ksp = ctx->ksp;
K = ctx->K;
G = ctx->G;
di = ctx->inv_diag_M;
s = ctx->s;
t = ctx->t;
X = ctx->X;
if (ctx->monitor_rhs_consistency) { BSSCRBSSCR_Lp_monitor_check_rhs_consistency(ksp,x,1); }
KSPSolve( ksp, x, t ); /* t <- GtG_inv x */
if (ctx->monitor_activated) { BSSCR_Lp_monitor(ksp,1); }
MatMult( G, t, s ); /* s <- G t */
VecPointwiseMult( s, s,di ); /* s <- s * di */
MatMultTranspose( K, s, X ); /* X <- K^T s */
VecPointwiseMult( X, X,di ); /* X <- X * di */
MatMultTranspose( G, X, t ); /* t <- Gt X */
if (ctx->monitor_rhs_consistency) { BSSCRBSSCR_Lp_monitor_check_rhs_consistency(ksp,t,2); }
KSPSolve( ksp, t, y ); /* y <- GtG_inv t */
if (ctx->monitor_activated) { BSSCR_Lp_monitor(ksp,2); }
PetscFunctionReturn(0);
}
开发者ID:AngelValverdePerez,项目名称:underworld2,代码行数:35,代码来源:pc_GtKG.c
示例5: PressurePoisson
PetscErrorCode PressurePoisson( UserContext *uc )
{
Mat A;
Vec b, px, py, p, u, v, ss, c;
KSP ksp;
PetscErrorCode ierr;
PetscFunctionBegin;
PetscLogEventBegin(EVENT_PressurePoisson,0,0,0,0);
PetscLogEventRegister(&EVENT_PressurePoisson,"PressurePoisson", 0);
/* Assemble and Solve for pressure */
ierr = PetscPrintf(PETSC_COMM_WORLD, "***********************\nPRESSURE\n"); CHKERRQ(ierr);
ierr = AssemblePressureMatrx( uc ); CHKERRQ(ierr);
ierr = AssemblePressureRHS( uc ); CHKERRQ(ierr);
ierr = SolvePressure( uc ); CHKERRQ(ierr);
ierr = WriteResult( uc, uc->p, "p_pressure" ); CHKERRQ(ierr);
/* Assemble and Solve for Velocity */
ierr = PetscPrintf(PETSC_COMM_WORLD, "***********************\nVELOCITY\n"); CHKERRQ(ierr);
ierr = AssembleVelocityRHS( uc ); CHKERRQ(ierr);
ierr = AssembleVelocityMatrix( uc ); CHKERRQ(ierr);
ierr = KSPSetOperators(uc->ksp,uc->A, uc->A, SAME_PRECONDITIONER); CHKERRQ(ierr);
ierr = KSPSolve(uc->ksp,uc->px,uc->u);CHKERRQ(ierr);
ierr = KSPSolve(uc->ksp,uc->py,uc->v);CHKERRQ(ierr);
ierr = WriteResult( uc, uc->u, "u_velocity" ); CHKERRQ(ierr);
ierr = WriteResult( uc, uc->v, "v_velocity" ); CHKERRQ(ierr);
ierr = ComputeShearStress(uc); CHKERRQ(ierr);
ierr = WriteResult( uc, uc->ss,"shear_stress" ); CHKERRQ(ierr);
ierr = ConservationTest(uc); CHKERRQ(ierr);
/* Output indexing */
PetscReal *idx;
VecGetArray(uc->b,&idx);
for( int i = 0; i < uc->numNodes; i++)
idx[i] = i;
VecRestoreArray(uc->b,&idx);
WriteResult(uc, uc->b, "indexes");
VecSet(uc->b, 0.);
ierr = VecDestroy(uc->c); CHKERRQ(ierr);
ierr = VecDestroy(uc->ss); CHKERRQ(ierr);
ierr = MatDestroy(uc->A); CHKERRQ(ierr);
ierr = KSPDestroy(uc->ksp); CHKERRQ(ierr);
ierr = VecDestroy(uc->b); CHKERRQ(ierr);
ierr = VecDestroy(uc->p); CHKERRQ(ierr);
ierr = VecDestroy(uc->px); CHKERRQ(ierr);
ierr = VecDestroy(uc->py); CHKERRQ(ierr);
ierr = VecDestroy(uc->u); CHKERRQ(ierr);
ierr = VecDestroy(uc->v); CHKERRQ(ierr);
PetscLogEventEnd(EVENT_PressurePoisson,0,0,0,0);
PetscFunctionReturn(0);
}
开发者ID:adrielb,项目名称:DCell,代码行数:56,代码来源:PressurePoisson.c
示例6: PCMGMCycle_Private
PetscErrorCode PCMGMCycle_Private(PC pc,PC_MG_Levels **mglevelsin,PCRichardsonConvergedReason *reason)
{
PC_MG *mg = (PC_MG*)pc->data;
PC_MG_Levels *mgc,*mglevels = *mglevelsin;
PetscErrorCode ierr;
PetscInt cycles = (mglevels->level == 1) ? 1 : (PetscInt) mglevels->cycles;
PC subpc;
PCFailedReason pcreason;
PetscFunctionBegin;
if (mglevels->eventsmoothsolve) {ierr = PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
ierr = KSPSolve(mglevels->smoothd,mglevels->b,mglevels->x);CHKERRQ(ierr); /* pre-smooth */
ierr = KSPGetPC(mglevels->smoothd,&subpc);CHKERRQ(ierr);
ierr = PCGetSetUpFailedReason(subpc,&pcreason);CHKERRQ(ierr);
if (pcreason) {
pc->failedreason = PC_SUBPC_ERROR;
}
if (mglevels->eventsmoothsolve) {ierr = PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
if (mglevels->level) { /* not the coarsest grid */
if (mglevels->eventresidual) {ierr = PetscLogEventBegin(mglevels->eventresidual,0,0,0,0);CHKERRQ(ierr);}
ierr = (*mglevels->residual)(mglevels->A,mglevels->b,mglevels->x,mglevels->r);CHKERRQ(ierr);
if (mglevels->eventresidual) {ierr = PetscLogEventEnd(mglevels->eventresidual,0,0,0,0);CHKERRQ(ierr);}
/* if on finest level and have convergence criteria set */
if (mglevels->level == mglevels->levels-1 && mg->ttol && reason) {
PetscReal rnorm;
ierr = VecNorm(mglevels->r,NORM_2,&rnorm);CHKERRQ(ierr);
if (rnorm <= mg->ttol) {
if (rnorm < mg->abstol) {
*reason = PCRICHARDSON_CONVERGED_ATOL;
ierr = PetscInfo2(pc,"Linear solver has converged. Residual norm %g is less than absolute tolerance %g\n",(double)rnorm,(double)mg->abstol);CHKERRQ(ierr);
} else {
*reason = PCRICHARDSON_CONVERGED_RTOL;
ierr = PetscInfo2(pc,"Linear solver has converged. Residual norm %g is less than relative tolerance times initial residual norm %g\n",(double)rnorm,(double)mg->ttol);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
}
mgc = *(mglevelsin - 1);
if (mglevels->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
ierr = MatRestrict(mglevels->restrct,mglevels->r,mgc->b);CHKERRQ(ierr);
if (mglevels->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
ierr = VecSet(mgc->x,0.0);CHKERRQ(ierr);
while (cycles--) {
ierr = PCMGMCycle_Private(pc,mglevelsin-1,reason);CHKERRQ(ierr);
}
if (mglevels->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
ierr = MatInterpolateAdd(mglevels->interpolate,mgc->x,mglevels->x,mglevels->x);CHKERRQ(ierr);
if (mglevels->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
if (mglevels->eventsmoothsolve) {ierr = PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
ierr = KSPSolve(mglevels->smoothu,mglevels->b,mglevels->x);CHKERRQ(ierr); /* post smooth */
if (mglevels->eventsmoothsolve) {ierr = PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
}
PetscFunctionReturn(0);
}
开发者ID:ziolai,项目名称:petsc,代码行数:56,代码来源:mg.c
示例7: PCApply_NN
static PetscErrorCode PCApply_NN(PC pc,Vec r,Vec z)
{
PC_IS *pcis = (PC_IS*)(pc->data);
PetscErrorCode ierr;
PetscScalar m_one = -1.0;
Vec w = pcis->vec1_global;
PetscFunctionBegin;
/*
Dirichlet solvers.
Solving $ B_I^{(i)}r_I^{(i)} $ at each processor.
Storing the local results at vec2_D
*/
ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecScatterEnd (pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = KSPSolve(pcis->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr);
/*
Computing $ r_B - \sum_j \tilde R_j^T A_{BI}^{(j)} (B_I^{(j)}r_I^{(j)}) $ .
Storing the result in the interface portion of the global vector w.
*/
ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr);
ierr = VecScale(pcis->vec1_B,m_one);CHKERRQ(ierr);
ierr = VecCopy(r,w);CHKERRQ(ierr);
ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,w,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,w,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
/*
Apply the interface preconditioner
*/
ierr = PCNNApplyInterfacePreconditioner(pc,w,z,pcis->work_N,pcis->vec1_B,pcis->vec2_B,pcis->vec3_B,pcis->vec1_D,
pcis->vec3_D,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr);
/*
Computing $ t_I^{(i)} = A_{IB}^{(i)} \tilde R_i z_B $
The result is stored in vec1_D.
*/
ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecScatterEnd (pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr);
/*
Dirichlet solvers.
Computing $ B_I^{(i)}t_I^{(i)} $ and sticking into the global vector the blocks
$ B_I^{(i)}r_I^{(i)} - B_I^{(i)}t_I^{(i)} $.
*/
ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
ierr = VecScatterEnd (pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
ierr = KSPSolve(pcis->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr);
ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr);
ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
ierr = VecScatterEnd (pcis->global_to_D,pcis->vec2_D,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:Kun-Qu,项目名称:petsc,代码行数:54,代码来源:nn.c
示例8: PCMGACycle_Private
PetscErrorCode PCMGACycle_Private(PC pc,PC_MG_Levels **mglevels)
{
PetscErrorCode ierr;
PetscInt i,l = mglevels[0]->levels;
PetscFunctionBegin;
/* compute RHS on each level */
for (i=l-1; i>0; i--) {
if (mglevels[i]->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
ierr = MatRestrict(mglevels[i]->restrct,mglevels[i]->b,mglevels[i-1]->b);CHKERRQ(ierr);
if (mglevels[i]->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
}
/* solve separately on each level */
for (i=0; i<l; i++) {
ierr = VecSet(mglevels[i]->x,0.0);CHKERRQ(ierr);
if (mglevels[i]->eventsmoothsolve) {ierr = PetscLogEventBegin(mglevels[i]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
ierr = KSPSolve(mglevels[i]->smoothd,mglevels[i]->b,mglevels[i]->x);CHKERRQ(ierr);
ierr = KSPCheckSolve(mglevels[i]->smoothd,pc,mglevels[i]->x);CHKERRQ(ierr);
if (mglevels[i]->eventsmoothsolve) {ierr = PetscLogEventEnd(mglevels[i]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);}
}
for (i=1; i<l; i++) {
if (mglevels[i]->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
ierr = MatInterpolateAdd(mglevels[i]->interpolate,mglevels[i-1]->x,mglevels[i]->x,mglevels[i]->x);CHKERRQ(ierr);
if (mglevels[i]->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);}
}
PetscFunctionReturn(0);
}
开发者ID:firedrakeproject,项目名称:petsc,代码行数:27,代码来源:smg.c
示例9: VecScale
bool PotentialSolve::Solve(Vec delta, Vec pot, double bias) {
PetscInt its;
KSPConvergedReason reason;
bool retval;
// Delta to -Delta
VecScale(delta, -1.0/bias);
// Actually solve
KSPSolve(solver, delta, pot);
KSPGetConvergedReason(solver,&reason);
if (reason<0) {
PetscPrintf(PETSC_COMM_WORLD,"Diverged : %d.\n",reason);
retval=false;
} else {
KSPGetIterationNumber(solver,&its);
PetscPrintf(PETSC_COMM_WORLD,"\nConvergence in %d iterations.\n",(int)its);
retval=true;
}
// Remove any values that might have crept in here
_dg1.ZeroPad(pot);
// Clean up
VecScale(delta, -1.0*bias);
return retval;
}
开发者ID:ajcuesta,项目名称:baorecon,代码行数:29,代码来源:PotentialSolve.cpp
示例10: swap
void SingleLongPipe::applyAdvection() {
swap(rhs, concentration);
PetscScalar* localRhs = rhs.getLocalArray();
PetscScalar alpha = flux*dt*T0/M_PI/(L*L*L);
const PetscScalar *localRadii = radii.getLocalArray();
PetscInt mStart, mEnd, jStart = 0;
MatGetOwnershipRange(M, &mStart, &mEnd);
if (rank == 0) {
int row = 0;
int column = 0;
PetscScalar val = alpha / (localRadii[0] * localRadii[0]);
PetscScalar val1 = 1 + val;
MatSetValues(M, 1, &row, 1, &column, &val1, INSERT_VALUES);
localRhs[0] += val*C0;
++mStart;
++jStart;
}
for (int i = mStart, j = jStart; i < mEnd; ++i, ++j) {
int row = i;
int column = i;
PetscScalar val = alpha / (localRadii[j] * localRadii[j]);
PetscScalar val1 = 1 + val;
MatSetValues(M, 1, &row, 1, &column, &val1, INSERT_VALUES);
val1 = -val;
--column;
MatSetValues(M, 1, &row, 1, &column, &val1, INSERT_VALUES);
}
MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
rhs.restoreVec(localRhs);
KSPSolve(ksp, rhs.getVec(), concentration.getVec());
}
开发者ID:whyzidane,项目名称:CO2_project,代码行数:32,代码来源:single_long_pipe_class.cpp
示例11: SFieldSolveFor
double SFieldSolveFor(SField sfv, double *Y, unsigned int yCount) {
mySField sf = static_cast<mySField>(sfv);
assert(yCount <= sf->maxN);
assert(Y);
assert(sf->running);
sf->Y = Y;
sf->curN = yCount;
// -------------- SOLVE
PetscErrorCode ierr;
PetscLogDouble tic,toc;
PetscTime(&tic);
int pt[sf->d];
ierr = MatZeroEntries(sf->J); CHKERRQ(ierr);
JacobianOnD(sf->J, sf->F, 0, pt, sf);
ierr = MatAssemblyBegin(sf->J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(sf->J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
PetscTime(&toc);
sf->timeAssembly += toc-tic;
PetscTime(&tic);
ierr = VecZeroEntries(sf->U); CHKERRQ(ierr);
ierr = KSPSetOperators(sf->ksp, sf->J, sf->J); CHKERRQ(ierr);
ierr = KSPSetUp(sf->ksp); CHKERRQ(ierr);
ierr = KSPSolve(sf->ksp,sf->F,sf->U); CHKERRQ(ierr);
PetscTime(&toc);
sf->timeSolver += toc-tic;
return Integrate(sf->U,pt,0,sf);
}
开发者ID:StochasticNumerics,项目名称:mimclib,代码行数:31,代码来源:matern.cpp
示例12: main
int main(int argc,char **argv)
{
PetscErrorCode ierr;
KSP ksp;
DM da,shell;
PetscInt levels;
PetscInitialize(&argc,&argv,(char*)0,help);
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-129,1,1,0,&da);CHKERRQ(ierr);
ierr = MyDMShellCreate(PETSC_COMM_WORLD,da,&shell);CHKERRQ(ierr);
/* these two lines are not needed but allow PCMG to automatically know how many multigrid levels the user wants */
ierr = DMGetRefineLevel(da,&levels);CHKERRQ(ierr);
ierr = DMSetRefineLevel(shell,levels);CHKERRQ(ierr);
ierr = KSPSetDM(ksp,shell);CHKERRQ(ierr);
ierr = KSPSetComputeRHS(ksp,ComputeRHS,NULL);CHKERRQ(ierr);
ierr = KSPSetComputeOperators(ksp,ComputeMatrix,NULL);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSolve(ksp,NULL,NULL);CHKERRQ(ierr);
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = DMDestroy(&shell);CHKERRQ(ierr);
ierr = DMDestroy(&da);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}
开发者ID:000Justin000,项目名称:ATPESC,代码行数:29,代码来源:ex65.c
示例13: MatAssemblyBegin
void PETSc::Solve(void)
{
//start_clock("Before Assemble matrix and vector");
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
ierr = VecAssemblyBegin(x);
ierr = VecAssemblyEnd(x);
ierr = VecAssemblyBegin(b);
ierr = VecAssemblyEnd(b);
//stop_clock("After Assembly matrix and vector");
KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
KSPSetType(ksp,KSPBCGSL);
KSPBCGSLSetEll(ksp,2);
//KSPGetPC(ksp, &pc);
//PCSetType(pc, PCJACOBI);
KSPSetFromOptions(ksp);
KSPSetUp(ksp);
//start_clock("Before KSPSolve");
KSPSolve(ksp,b,x);
//stop_clock("After KSPSolve");
}
开发者ID:tonyguo1,项目名称:LSGFD,代码行数:29,代码来源:solver_petsc.cpp
示例14: SolvePressure
PetscErrorCode SolvePressure(UserContext *uc)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscLogEventBegin(EVENT_SolvePressure,0,0,0,0);
ierr = KSPCreate(PETSC_COMM_WORLD, &uc->ksp); CHKERRQ(ierr);
ierr = KSPSetType(uc->ksp,KSPCG); CHKERRQ(ierr);
ierr = KSPSetOperators(uc->ksp,uc->A,uc->A,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
ierr = KSPSetFromOptions(uc->ksp); CHKERRQ(ierr);
ierr = KSPSetTolerances(uc->ksp, 0.000001, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT); CHKERRQ(ierr);
ierr = KSPSetType( uc->ksp, KSPPREONLY); CHKERRQ(ierr);
PC pc;
ierr = KSPGetPC(uc->ksp, &pc); CHKERRQ(ierr);
ierr = PCSetType(pc, PCLU); CHKERRQ(ierr);
ierr = KSPSolve(uc->ksp,uc->b,uc->p);CHKERRQ(ierr);
KSPConvergedReason reason;
KSPGetConvergedReason(uc->ksp,&reason);
PetscPrintf(PETSC_COMM_WORLD,"Pressure KSPConvergedReason: %D\n", reason);
// if( reason < 0 ) SETERRQ(PETSC_ERR_CONV_FAILED, "Exiting: Failed to converge\n");
PetscLogEventEnd(EVENT_SolvePressure,0,0,0,0);
PetscFunctionReturn(0);
}
开发者ID:adrielb,项目名称:DCell,代码行数:28,代码来源:PressurePoisson.c
示例15: main
int main(int argc, char **args) {
PetscErrorCode ierr;
ierr = DCellInit(); CHKERRQ(ierr);
PetscReal dx = 1;
iCoor size = {1625,1145,0};
// iCoor size = {253,341,0};
int fd;
Coor dh = {dx,dx,dx};
iCoor pos = {0,0,0};
Grid chip;
ierr = GridCreate(dh,pos,size,1,&chip); CHKERRQ(ierr);
ierr = PetscInfo(0,"Reading image file\n"); CHKERRQ(ierr);
ierr = PetscBinaryOpen("/scratch/n/BL Big",FILE_MODE_READ,&fd); CHKERRQ(ierr);
ierr = PetscBinaryRead(fd,chip->v1,size.x*size.y,PETSC_DOUBLE); CHKERRQ(ierr);
ierr = PetscBinaryClose(fd); CHKERRQ(ierr);
ierr = PetscInfo(0,"Writing image file\n"); CHKERRQ(ierr);
ierr = GridWrite(chip,0); CHKERRQ(ierr);
FluidField fluid;
ierr = FluidFieldCreate(PETSC_COMM_WORLD, &fluid); CHKERRQ(ierr);
ierr = FluidFieldSetDims(fluid,size); CHKERRQ(ierr);
ierr = FluidFieldSetDx(fluid,dx); CHKERRQ(ierr);
ierr = FluidFieldSetMask(fluid, chip); CHKERRQ(ierr);
ierr = FluidFieldSetup(fluid); CHKERRQ(ierr);
ierr = SetPressureBC(fluid); CHKERRQ(ierr);
ierr = KSPSolve(fluid->ksp,fluid->rhs,fluid->vel); CHKERRQ(ierr);
ierr = FluidFieldWrite( fluid,0); CHKERRQ(ierr);
ierr = FluidFieldDestroy(fluid); CHKERRQ(ierr);
ierr = GridDestroy(chip); CHKERRQ(ierr);
ierr = DCellFinalize(); CHKERRQ(ierr);
return 0;
}
开发者ID:adrielb,项目名称:DCell,代码行数:34,代码来源:Microfluidics.c
示例16: main
int main(int argc,char **argv)
{
KSP ksp;
DM da;
UserContext user;
PetscInt bc;
PetscErrorCode ierr;
PetscInitialize(&argc,&argv,(char *)0,help);
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-11,-11,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&da);CHKERRQ(ierr);
ierr = KSPSetDM(ksp,(DM)da);
ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr);
user.uu = 1.0;
user.tt = 1.0;
bc = (PetscInt)NEUMANN; // Use Neumann Boundary Conditions
user.bcType = (BCType)bc;
ierr = KSPSetComputeRHS(ksp,ComputeRHS,&user);CHKERRQ(ierr);
ierr = KSPSetComputeOperators(ksp,ComputeJacobian,&user);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSolve(ksp,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
ierr = DMDestroy(&da);CHKERRQ(ierr);
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = PetscFinalize();CHKERRQ(ierr);
return 0;
}
开发者ID:erdc-cm,项目名称:petsc-dev,代码行数:30,代码来源:ex50.c
示例17: START_TEST
END_TEST
START_TEST( Generate2DLapacian_test )
{
PetscErrorCode ierr;
PetscInt d = 50, len = d*d;
Mat m;
PetscViewer view;
Generate2DLapacian( d, d, &m);
PetscViewerASCIIOpen(PETSC_COMM_SELF, "mat.dat", &view);
MatView(m, view);
KSP ksp;
PC pc;
ierr = KSPCreate(PETSC_COMM_SELF, &ksp); CHKERRQ(ierr);
KSPSetType(ksp, KSPPREONLY);
KSPGetPC(ksp, &pc);
PCSetType(pc, PCCHOLESKY);
KSPSetOperators(ksp, m, m, SAME_PRECONDITIONER);
Vec v, s;
VecCreateSeq(PETSC_COMM_SELF,len,&v);
VecDuplicate(v, &s);
VecSet(v, 1);
ierr = KSPSolve(ksp, v, s); CHKERRQ(ierr);
}
开发者ID:adrielb,项目名称:DCell,代码行数:27,代码来源:check.c
示例18: _check_if_constant_nullsp_present
Bool _check_if_constant_nullsp_present( Stokes_SLE_UzawaSolver* self, Mat K, Mat G, Mat M, Vec t1, Vec ustar, Vec r, Vec l, KSP ksp )
{
PetscInt N;
PetscScalar sum;
PetscReal nrm;
Bool nullsp_present;
VecGetSize(l,&N);
sum = 1.0/N;
VecSet(l,sum);
/* [S] {l} = {r} */
MatMult( G,l, t1 );
KSPSolve( ksp, t1, ustar );
MatMultTranspose( G, ustar, r );
if ( M ) {
VecScale( r, -1.0 );
MatMultAdd( M,l, r, r );
VecScale( r, -1.0 );
}
VecNorm(r,NORM_2,&nrm);
if (nrm < 1.e-7) {
Journal_PrintfL( self->info, 1, "Constant null space detected, " );
nullsp_present = True;
}
else {
Journal_PrintfL( self->info, 1, "Constant null space not present, " );
nullsp_present = False;
}
Journal_PrintfL( self->info, 1, "|| [S]{1} || = %G\n", nrm );
return nullsp_present;
}
开发者ID:OlympusMonds,项目名称:EarthByte_Underworld,代码行数:35,代码来源:Stokes_SLE_UzawaSolver.c
示例19: _Stokes_SLE_UzawaSolver_GetRhs
PetscErrorCode _Stokes_SLE_UzawaSolver_GetRhs( void *stokesSLE, void *solver, Vec rhs )
{
Stokes_SLE_UzawaSolver *self = (Stokes_SLE_UzawaSolver*)solver;
Stokes_SLE *sle = (Stokes_SLE*)stokesSLE;
/* stg linear algebra */
KSP A11_solver = self->velSolver;
Mat A12 = sle->gStiffMat->matrix;
Vec b1 = sle->fForceVec->vector;
Vec b2 = sle->hForceVec->vector;
Vec u_star = self->vStarVec;
/* petsc variables */
KSP ksp_A11;
/* check operations will be valid */
if (sle->dStiffMat!=NULL) { Stg_SETERRQ( PETSC_ERR_SUP, "A21 must be NULL" ); }
if (A11_solver==NULL){ Stg_SETERRQ( PETSC_ERR_ARG_NULL, "vel_solver is NULL" ); }
if (A12==NULL){ Stg_SETERRQ( PETSC_ERR_ARG_NULL, "A12 is NULL" ); }
if (b1==NULL){ Stg_SETERRQ( PETSC_ERR_ARG_NULL, "b1 is NULL" ); }
if (b2==NULL){ Stg_SETERRQ( PETSC_ERR_ARG_NULL, "b2 is NULL" ); }
if (u_star==NULL){ Stg_SETERRQ( PETSC_ERR_ARG_NULL, "u* is NULL" ); }
/* Extract petsc objects */
ksp_A11 = A11_solver;
/* compute rhs = A12^T A11^{-1} b1 - b2 */
KSPSolve( ksp_A11, b1, u_star ); /* u* = A11^{-1} b1 */
MatMultTranspose( A12, u_star, rhs ); /* b2 = A12^T u* */
VecAXPY( rhs, -1.0, b2 ); /* rhs <- rhs - b2 */
}
开发者ID:OlympusMonds,项目名称:EarthByte_Underworld,代码行数:30,代码来源:Stokes_SLE_UzawaSolver.c
示例20: PCApply_Redistribute
static PetscErrorCode PCApply_Redistribute(PC pc,Vec b,Vec x)
{
PC_Redistribute *red = (PC_Redistribute*)pc->data;
PetscErrorCode ierr;
PetscInt dcnt = red->dcnt,i;
const PetscInt *drows = red->drows;
PetscScalar *xwork;
const PetscScalar *bwork,*diag = red->diag;
PetscFunctionBegin;
if (!red->work) {
ierr = VecDuplicate(b,&red->work);CHKERRQ(ierr);
}
/* compute the rows of solution that have diagonal entries only */
ierr = VecSet(x,0.0);CHKERRQ(ierr); /* x = diag(A)^{-1} b */
ierr = VecGetArray(x,&xwork);CHKERRQ(ierr);
ierr = VecGetArrayRead(b,&bwork);CHKERRQ(ierr);
for (i=0; i<dcnt; i++) {
xwork[drows[i]] = diag[i]*bwork[drows[i]];
}
ierr = PetscLogFlops(dcnt);CHKERRQ(ierr);
ierr = VecRestoreArray(red->work,&xwork);CHKERRQ(ierr);
ierr = VecRestoreArrayRead(b,&bwork);CHKERRQ(ierr);
/* update the right hand side for the reduced system with diagonal rows (and corresponding columns) removed */
ierr = MatMult(pc->pmat,x,red->work);CHKERRQ(ierr);
ierr = VecAYPX(red->work,-1.0,b);CHKERRQ(ierr); /* red->work = b - A x */
ierr = VecScatterBegin(red->scatter,red->work,red->b,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = VecScatterEnd(red->scatter,red->work,red->b,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
ierr = KSPSolve(red->ksp,red->b,red->x);CHKERRQ(ierr);
ierr = VecScatterBegin(red->scatter,red->x,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
ierr = VecScatterEnd(red->scatter,red->x,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
开发者ID:erdc-cm,项目名称:petsc-dev,代码行数:34,代码来源:redistribute.c
注:本文中的KSPSolve函数示例由纯净天空整理自Github/MSDocs等源码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。 |
请发表评论