Skip to content

Commit

Permalink
In progress of changing structure SlaybaughLab#10
Browse files Browse the repository at this point in the history
  • Loading branch information
Weixiong Zheng committed Aug 21, 2017
1 parent 2bfab7c commit f0e4705
Show file tree
Hide file tree
Showing 20 changed files with 458 additions and 223 deletions.
11 changes: 5 additions & 6 deletions src/common/bart_driver.cc
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,7 @@ void BartDriver<dim>::initialize_dealii_objects ()
MPI_COMM_WORLD,
relevant_dofs);

itr_cls.initialize_system_matrices_vectors (dsp);
itr_cls.initialize_system_matrices_vectors (dsp, local_dofs);
}

template <int dim>
Expand All @@ -129,8 +129,7 @@ void BartDriver<dim>::output_results () const
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);

if (zeroth_proc())
itr_cls.get_keff (keff);
itr_cls.get_keff (keff);
itr_cls.get_flux_this_proc (sflx_proc);

data_out.add_data_vector (keff, "keff");
Expand All @@ -153,7 +152,7 @@ void BartDriver<dim>::output_results () const
std::ofstream output ((filename + ".vtu").c_str ());
data_out.write_vtu (output);

if (zeroth_proc())
if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0)
{
std::vector<std::string> filenames;
for (unsigned int i=0;
Expand All @@ -171,11 +170,11 @@ void BartDriver<dim>::output_results () const
template <int dim>
void BartDriver<dim>::run ()
{
radio ("making grid");
msh_ptr->make_grid (triangulation);
setup_system ();
report_system ();
itr_cls.do_iterations (msh_ptr, aqd_ptr);
itr_cls.solve_problems (msh_ptr, aqd_ptr, sflx_proc);
itr_cls.get_keff (keff);
output_results ();
}

Expand Down
1 change: 0 additions & 1 deletion src/common/bart_driver.h
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,6 @@ class BartDriver
std_cxx11::shared_ptr<MeshGenerator<dim> > msh_ptr;
std_cxx11::shared_ptr<MaterialProperties> mat_ptr;
std_cxx11::shared_ptr<AQBase<dim> > aqd_ptr;
std_cxx11::shared_ptr<PreconditionerSolver> sol_ptr;

std::string transport_model_name;
std::string ho_linear_solver_name;
Expand Down
11 changes: 6 additions & 5 deletions src/common/bart_tools.cc
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,12 @@ FE_Poly<TensorProductPolynomials<dim>,dim,dim>* build_finite_element
return new FE_Q<dim> (p_order);
}

template <int dim>
std_cxx11::shared_ptr<PreconditionerSolver> build_linalg (ParameterHandler &prm)
{
return std_cxx11::shared_ptr<PreconditionerSolver> (new PreconditionerSolver(prm));
}

template <int dim>
std_cxx11::shared_ptr<MeshGenerator<dim> > build_mesh (ParameterHandler &prm)
{
Expand Down Expand Up @@ -66,11 +72,6 @@ build_aq_model (ParameterHandler &prm)
return aq_class;
}

bool zeroth_proc ()
{
return Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0
}

void radio (std::string str)
{
pout << str << std::endl;
Expand Down
124 changes: 60 additions & 64 deletions src/common/preconditioner_solver.cc
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ void PreconditionerSolver::initialize_ho_preconditioners
data.symmetric_operator = false;
else
data.symmetric_operator = true;
pre_ho_amg[i]->initialize(*(ho_syses)[i], data);
pre_ho_amg[i]->initialize(*ho_syses[i], data);
}
}
else if (ho_preconditioner_name=="bjacobi")
Expand All @@ -66,7 +66,7 @@ void PreconditionerSolver::initialize_ho_preconditioners
{
pre_ho_bjacobi[i] = std_cxx11::shared_ptr<PETScWrappers::PreconditionBlockJacobi>
(new PETScWrappers::PreconditionBlockJacobi);
pre_ho_bjacobi[i]->initialize(*(ho_syses)[i]);
pre_ho_bjacobi[i]->initialize(*ho_syses[i]);
}
}
else if (ho_preconditioner_name=="jacobi")
Expand All @@ -76,7 +76,7 @@ void PreconditionerSolver::initialize_ho_preconditioners
{
pre_ho_jacobi[i] = std_cxx11::shared_ptr<PETScWrappers::PreconditionJacobi>
(new PETScWrappers::PreconditionJacobi);
pre_ho_jacobi[i]->initialize(*(ho_syses)[i]);
pre_ho_jacobi[i]->initialize(*ho_syses[i]);
}
}
else if (ho_preconditioner_name=="bssor")
Expand All @@ -87,7 +87,7 @@ void PreconditionerSolver::initialize_ho_preconditioners
pre_ho_eisenstat[i] = std_cxx11::shared_ptr<PETScWrappers::PreconditionEisenstat>
(new PETScWrappers::PreconditionEisenstat);
PETScWrappers::PreconditionEisenstat::AdditionalData data(ho_ssor_omega);
pre_ho_eisenstat[i]->initialize(*(ho_syses)[i], data);
pre_ho_eisenstat[i]->initialize(*ho_syses[i], data);
}
}
else if (ho_preconditioner_name=="parasails")
Expand All @@ -101,12 +101,12 @@ void PreconditionerSolver::initialize_ho_preconditioners
(transport_model_name=="ep" && have_reflective_bc))
{
PETScWrappers::PreconditionParaSails::AdditionalData data (2);
pre_ho_parasails[i]->initialize(*(ho_syses)[i], data);
pre_ho_parasails[i]->initialize(*ho_syses[i], data);
}
else
{
PETScWrappers::PreconditionParaSails::AdditionalData data (1);
pre_ho_parasails[i]->initialize(*(ho_syses)[i], data);
pre_ho_parasails[i]->initialize(*ho_syses[i], data);
}
}
}
Expand All @@ -124,89 +124,85 @@ void PreconditionerSolver::initialize_ho_preconditioners
}

void PreconditionerSolver::ho_solve
(std::vector<PETScWrappers::MPI::SparseMatrix*> &ho_syses,
std::vector<PETScWrappers::MPI::Vector*> &ho_psis,
std::vector<PETScWrappers::MPI::Vector*> &ho_rhses,
(PETScWrappers::MPI::SparseMatrix &ho_sys,
PETScWrappers::MPI::Vector &ho_psi,
PETScWrappers::MPI::Vector &ho_rhs,
unsigned int i/*component number*/)
{
AssertThrow (n_total_ho_vars==ho_syses.size(),
ExcMessage("num of HO system matrices should be equal to total variable number"));
AssertThrow (n_total_ho_vars==ho_rhses.size(),
ExcMessage("num of HO system rhs should be equal to total variable number"));
if (ho_linear_solver_name=="cg")
{
PETScWrappers::SolverCG
solver (*ho_cn[i], MPI_COMM_WORLD);
if (ho_preconditioner_name=="amg")
solver.solve (*ho_syses[i],
*ho_psis[i],
*ho_rhses[i],
*(pre_ho_amg)[i]);
solver.solve (ho_sys,
ho_psi,
ho_rhs,
*pre_ho_amg[i]);
else if (ho_preconditioner_name=="jacobi")
solver.solve (*ho_syses[i],
*ho_psis[i],
*ho_rhses[i],
*(pre_ho_jacobi)[i]);
solver.solve (ho_sys,
ho_psi,
ho_rhs,
*pre_ho_jacobi[i]);
else if (ho_preconditioner_name=="bssor")
solver.solve (*ho_syses[i],
*ho_psis[i],
*ho_rhses[i],
*(pre_ho_eisenstat)[i]);
solver.solve (ho_sys,
ho_psi,
ho_rhs,
*pre_ho_eisenstat[i]);
else if (ho_preconditioner_name=="parasails")
solver.solve (*ho_syses[i],
*ho_psis[i],
*ho_rhses[i],
*(pre_ho_parasails)[i]);
solver.solve (ho_sys,
ho_psi,
ho_rhs,
*pre_ho_parasails[i]);
}
else if (ho_linear_solver_name=="bicgstab")
{
PETScWrappers::SolverBicgstab
solver (*ho_cn[i], MPI_COMM_WORLD);
if (ho_preconditioner_name=="amg")
solver.solve (*ho_syses[i],
*ho_psis[i],
*ho_rhses[i],
*(pre_ho_amg)[i]);
solver.solve (ho_sys,
ho_psi,
ho_rhs,
*pre_ho_amg[i]);
else if (ho_preconditioner_name=="jacobi")
solver.solve (*ho_syses[i],
*ho_psis[i],
*ho_rhses[i],
*(pre_ho_jacobi)[i]);
solver.solve (ho_sys,
ho_psi,
ho_rhs,
*pre_ho_jacobi[i]);
else if (ho_preconditioner_name=="bssor")
solver.solve (*ho_syses[i],
*ho_psis[i],
*ho_rhses[i],
*(pre_ho_eisenstat)[i]);
solver.solve (ho_sys,
ho_psi,
ho_rhs,
*pre_ho_eisenstat[i]);
else if (ho_preconditioner_name=="parasails")
solver.solve (*ho_syses[i],
*ho_psis[i],
*ho_rhses[i],
*(pre_ho_parasails)[i]);
solver.solve (ho_sys,
ho_psi,
ho_rhs,
*pre_ho_parasails[i]);
}
else if (ho_linear_solver_name=="gmres")
{
PETScWrappers::SolverGMRES
solver (*ho_cn[i], MPI_COMM_WORLD);
if (ho_preconditioner_name=="amg")
solver.solve (*ho_syses[i],
*ho_psis[i],
*ho_rhses[i],
*(pre_ho_amg)[i]);
solver.solve (ho_sys,
ho_psi,
ho_rhs,
*pre_ho_amg[i]);
else if (ho_preconditioner_name=="jacobi")
solver.solve (*ho_syses[i],
*ho_psis[i],
*ho_rhses[i],
*(pre_ho_jacobi)[i]);
solver.solve (ho_sys,
ho_psi,
ho_rhs,
*pre_ho_jacobi[i]);
else if (ho_preconditioner_name=="bssor")
solver.solve (*ho_syses[i],
*ho_psis[i],
*ho_rhses[i],
*(pre_ho_eisenstat)[i]);
solver.solve (ho_sys,
ho_psi,
ho_rhs,
*pre_ho_eisenstat[i]);
else if (ho_preconditioner_name=="parasails")
solver.solve (*ho_syses[i],
*ho_psis[i],
*ho_rhses[i],
*(pre_ho_parasails)[i]);
solver.solve (ho_sys,
ho_psi,
ho_rhs,
*pre_ho_parasails[i]);
}
else// if (linear_solver_name=="direct")
{
Expand All @@ -221,9 +217,9 @@ void PreconditionerSolver::ho_solve
ho_direct[i]->set_symmetric_mode (true);
ho_direct_init[i] = true;
}
ho_direct[i]->solve (*ho_syses[i],
*ho_psis[i],
*ho_rhses[i]);
ho_direct[i]->solve (ho_sys,
ho_psi,
ho_rhs);
}
// the ho_linear_iters are for reporting linear solver status, test purpose only
if (ho_linear_solver_name!="direct")
Expand Down
62 changes: 62 additions & 0 deletions src/iteration/eigen_base.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,62 @@
#include "eigen_base.h"

template <int dim>
EigenBase<dim>::EigenBase () : IterationBase<dim> (),
err_k_tol(1.0e-6),
err_phi_tol(1.0e-7),
err_phi_eigen_tol(1.0e-5),
keff(1.0)
{
}

template <int dim>
EigenBase<dim>::~EigenBase ()
{
}

template <int dim>
EigenBase<dim>::do_iterations ()
{
this->initialize_equations ();
eigen_iterations ();
}

template <int dim>
void EigenBase<dim>::generate_system_matrices
(std::vector<PETScWrappers::MPI::SparseMatrix*> sys_mats)
{
// EquationBase will be instantiated in InGroupBase
mgs_ptr->generate_system_matrices (sys_mats);
}

template <int dim>
void EigenBase<dim>::initialize_fiss_process ()
{
for (unsigned int g=0; g<n_group; ++g)
this->sflx_proc[g] = 1.0;

fission_source = this->trm_ptr->estimate_fiss_source (this->sflx_proc);
}

template <int dim>
double EigenBase<dim>::estimate_k (double &fiss_source,
double &fiss_source_prev,
double &k_prev)
{
return k_prev * fiss_source / fiss_source_prev;
}

template <int dim>
double EigenBase<dim>::estimate_k_err (double &k, double &k_prev)
{
return std::fabs (k - k_prev)/k;
}

template <int dim>
void EigenBase<dim>::eigen_iteration (double &keff,
std_cxx11::shared_ptr<MGSolver<dim> > mgs_ptr)
{
}

template class EigenBase<2>;
template class EigenBase<3>;
25 changes: 25 additions & 0 deletions src/iteration/eigen_base.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,25 @@
#ifndef __eigen_base_h__
#define __eigen_base_h__

template <int dim>
class EigenBase : public IterationBase
{
public:
EigenBase ();
virtual ~EigenBase ();

virtual void do_iterations ();
virtual void eigen_iterations ();

double get_keff ();

protected:
double estimate_fission_source (std::vector<Vector<double> > &sflx_proc);
double estimate_k (double &fiss_src, double &fiss_src_prev_gen, double &k_prev);
double estimate_k_err (double &k, double &k_prev);

double keff;
double keff_prev;
}

#endif //__eigen_base_h__
23 changes: 23 additions & 0 deletions src/iteration/in_group_base.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
#include "in_group_base.h"
#include "../common/preconditioner_solver.h"

template <int dim>
InGroupBase<dim>::InGroupBase ()
:
IterationBase<dim> ()
{
sol_ptr = build_linalg (prm);
}

template <int dim>
InGroupBase<dim>::~InGroupBase ()
{
}

template <int dim>
InGroupBase<dim>::solve_in_group ()
{
}

template class InGroupBase<2>;
template class InGroupBase<3>;
Loading

0 comments on commit f0e4705

Please sign in to comment.