C++ KSPSolve函数代码示例

本文整理汇总了C++中KSPSolve函数的典型用法代码示例。如果您正苦于以下问题：C++ KSPSolve函数的具体用法？C++ KSPSolve怎么用？C++ KSPSolve使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。

在下文中一共展示了KSPSolve函数的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的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

注：本文中的KSPSolve函数示例由纯净天空整理自Github/MSDocs等开源代码及文档管理平台，相关代码片段筛选自各路编程大神贡献的开源项目，源码版权归原作者所有，传播和使用请参考对应项目的License；未经允许，请勿转载。