packages/Alea/mc_solvers/LinearSystemFactory.cc
//---------------------------------------------------------------------------//
// Build matrix and RHS for diffusion problem
//---------------------------------------------------------------------------//
Build matrix and RHS for diffusion problem
Build matrix and RHS for diffusion problem
void LinearSystemFactory::buildConvectionDiffusionSystem(
    Teuchos::RCP<Teuchos::ParameterList> pl,
    Teuchos::RCP<CRS_MATRIX>                 &A,
    Teuchos::RCP<MV>                     &b )
{
    GO N = pl->get("matrix_size",25);
    SCALAR nu = pl->get("viscosity",1.0);
    Teuchos::RCP<CRS_MATRIX> Ad = buildLaplacianMatrix(N,pl);
    Teuchos::RCP<CRS_MATRIX> Ac = buildConvectionMatrix(N,pl);
    A = Teuchos::rcp_dynamic_cast<CRS_MATRIX>(
        Ad->add(1.0,*Ac,nu,Ad->getDomainMap(),Ad->getRangeMap(),pl),true);
    CHECK( A != Teuchos::null );
    b = Teuchos::rcp( new MV(A->getDomainMap(),1) );
    const double pi = 4.0*std::atan(1.0);
    Teuchos::ArrayRCP<double> b_data = b->getDataNonConst(0);
    LO Nlocal = b_data.size();
    for( LO i=0; i<Nlocal; ++i )
    {
        double gid = static_cast<double>(b->getMap()->getGlobalElement(i));
        double val = std::sin(2*pi*gid/static_cast<double>(N-1));
        b->replaceLocalValue(i,0,std::sin(val));
    }
}
void LinearSystemFactory :: buildConvectionDiffusionSystem ( Teuchos :: RCP < Teuchos :: ParameterList > pl , Teuchos :: RCP < CRS_MATRIX > & A , Teuchos :: RCP < MV > & b ) { GO N = pl -> get ( " " , 25 ) ; SCALAR nu = pl -> get ( " " , 1.0 ) ; Teuchos :: RCP < CRS_MATRIX > Ad = buildLaplacianMatrix ( N , pl ) ; Teuchos :: RCP < CRS_MATRIX > Ac = buildConvectionMatrix ( N , pl ) ; A = Teuchos :: rcp_dynamic_cast < CRS_MATRIX > ( Ad -> add ( 1.0 , * Ac , nu , Ad -> getDomainMap ( ) , Ad -> getRangeMap ( ) , pl ) , true ) ; CHECK ( A != Teuchos :: null ) ; b = Teuchos :: rcp ( new MV ( A -> getDomainMap ( ) , 1 ) ) ; const double pi = 4.0 * std :: atan ( 1.0 ) ; Teuchos :: ArrayRCP < double > b_data = b -> getDataNonConst ( 0 ) ; LO Nlocal = b_data . size ( ) ; for ( LO i = 0 ; i < Nlocal ; ++ i ) { double gid = static_cast < double > ( b -> getMap ( ) -> getGlobalElement ( i ) ) ; double val = std :: sin ( 2 * pi * gid / static_cast < double > ( N - 1 ) ) ; b -> replaceLocalValue ( i , 0 , std :: sin ( val ) ) ; } }
Build matrix and RHS for diffusion problem
Build matrix and RHS for diffusion problem
// Sinusoidal source with homogeneous Dirichlet boundaries