Skip to content

Commit

Permalink
Change TransportBase to EquationBase SlaybaughLab#10
Browse files Browse the repository at this point in the history
  • Loading branch information
Weixiong Zheng committed Aug 5, 2017
1 parent 0046d93 commit 1bdf3d8
Show file tree
Hide file tree
Showing 9 changed files with 103 additions and 246 deletions.
40 changes: 17 additions & 23 deletions src/common/bart_driver.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@

#include <algorithm>

#include "transport_base.h"
#include "equation_base.h"
#include "../aqdata/aq_base.h"
#include "../aqdata/aq_lsgc.h"

Expand All @@ -23,10 +23,8 @@ triangulation (MPI_COMM_WORLD,
(Triangulation<dim>::smoothing_on_refinement |
Triangulation<dim>::smoothing_on_coarsening)),
dof_handler (triangulation),
prm(prm),
transport_model_name(prm.get("transport model")),
aq_name(prm.get("angular quadrature name")),
discretization(prm.get("spatial discretization")),
n_group(prm.get_integer("number of groups")),
n_azi(prm.get_integer("angular quadrature order")),
is_eigen_problem(prm.get_bool("do eigenvalue calculations")),
Expand All @@ -47,6 +45,7 @@ pcout(std::cout,
n_dir = aqd_ptr->get_n_dir ();
msh_ptr = build_mesh (prm);
itr_ptr = build_transport_iteration (prm);
fe = build_finite_element (prm);
}

template <int dim>
Expand Down Expand Up @@ -130,22 +129,22 @@ void BartDriver<dim>::initialize_system_matrices_vectors ()
if (do_nda)
{
vec_lo_sys.push_back (new LA::MPI::SparseMatrix);
vec_lo_rhs.push_back (new LA::MPI::Vector);
//vec_lo_rhs.push_back (new LA::MPI::Vector);
vec_lo_sflx.push_back (new LA::MPI::Vector);
vec_lo_sflx_old.push_back (new LA::MPI::Vector);
vec_lo_fixed_rhs.push_back (new LA::MPI::Vector);
//vec_lo_sflx_old.push_back (new LA::MPI::Vector);
//vec_lo_fixed_rhs.push_back (new LA::MPI::Vector);
}

vec_ho_sflx.push_back (new LA::MPI::Vector);
vec_ho_sflx_prev_gen.push_back (new LA::MPI::Vector);
vec_ho_sflx_old.push_back (new LA::MPI::Vector);
//vec_ho_sflx_prev_gen.push_back (new LA::MPI::Vector);
//vec_ho_sflx_old.push_back (new LA::MPI::Vector);

for (unsigned int i_dir=0; i_dir<n_dir; ++i_dir)
{
vec_ho_sys.push_back (new LA::MPI::SparseMatrix);
vec_aflx.push_back (new LA::MPI::Vector);
vec_ho_rhs.push_back (new LA::MPI::Vector);
vec_ho_fixed_rhs.push_back (new LA::MPI::Vector);
//vec_aflx.push_back (new LA::MPI::Vector);
//vec_ho_rhs.push_back (new LA::MPI::Vector);
//vec_ho_fixed_rhs.push_back (new LA::MPI::Vector);
}
}

Expand All @@ -157,14 +156,14 @@ void BartDriver<dim>::initialize_system_matrices_vectors ()
local_dofs,
dsp,
MPI_COMM_WORLD);
vec_lo_rhs[g]->reinit (local_dofs, MPI_COMM_WORLD);
vec_lo_fixed_rhs[g]->reinit (local_dofs, MPI_COMM_WORLD);
//vec_lo_rhs[g]->reinit (local_dofs, MPI_COMM_WORLD);
//vec_lo_fixed_rhs[g]->reinit (local_dofs, MPI_COMM_WORLD);
vec_lo_sflx[g]->reinit (local_dofs, MPI_COMM_WORLD);
vec_lo_sflx_old[g]->reinit (local_dofs, MPI_COMM_WORLD);
//vec_lo_sflx_old[g]->reinit (local_dofs, MPI_COMM_WORLD);
}

vec_ho_sflx[g]->reinit (local_dofs, MPI_COMM_WORLD);
vec_ho_sflx_old[g]->reinit (local_dofs, MPI_COMM_WORLD);
//vec_ho_sflx_old[g]->reinit (local_dofs, MPI_COMM_WORLD);
}

for (unsigned int k=0; k<n_total_ho_vars; ++k)
Expand All @@ -173,20 +172,15 @@ void BartDriver<dim>::initialize_system_matrices_vectors ()
local_dofs,
dsp,
MPI_COMM_WORLD);
vec_aflx[k]->reinit(local_dofs, MPI_COMM_WORLD);
vec_ho_rhs[k]->reinit (local_dofs, MPI_COMM_WORLD);
vec_ho_fixed_rhs[k]->reinit (local_dofs, MPI_COMM_WORLD);
//vec_aflx[k]->reinit(local_dofs, MPI_COMM_WORLD);
//vec_ho_rhs[k]->reinit (local_dofs, MPI_COMM_WORLD);
//vec_ho_fixed_rhs[k]->reinit (local_dofs, MPI_COMM_WORLD);
}
}

template <int dim>
void BartDriver<dim>::initialize_dealii_objects ()
{
if (discretization=="dfem")
fe = (new FE_DGQ<dim> (p_order));
else
fe = (new FE_Q<dim> (p_order));

dof_handler.distribute_dofs (*fe);

local_dofs = dof_handler.locally_owned_dofs ();
Expand Down
56 changes: 3 additions & 53 deletions src/common/bart_driver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,63 +43,13 @@
using namespace dealii;

template <int dim>
class TransportBase
class BartDriver
{
public:
TransportBase (ParameterHandler &prm);// : ProblemDefinition<dim> (prm){}
virtual ~TransportBase ();
BartDriver (ParameterHandler &prm);// : ProblemDefinition<dim> (prm){}
virtual ~BartDriver ();

void run ();

virtual void pre_assemble_cell_matrices
(const std_cxx11::shared_ptr<FEValues<dim> > fv,
typename DoFHandler<dim>::active_cell_iterator &cell,
std::vector<std::vector<FullMatrix<double> > > &streaming_at_qp,
std::vector<FullMatrix<double> > &collision_at_qp);

virtual void integrate_cell_bilinear_form
(const std_cxx11::shared_ptr<FEValues<dim> > fv,
typename DoFHandler<dim>::active_cell_iterator &cell,
FullMatrix<double> &cell_matrix,
unsigned int &i_dir,
unsigned int &g,
std::vector<std::vector<FullMatrix<double> > > &streaming_at_qp,
std::vector<FullMatrix<double> > &collision_at_qp);

virtual void integrate_boundary_bilinear_form
(const std_cxx11::shared_ptr<FEFaceValues<dim> > fvf,
typename DoFHandler<dim>::active_cell_iterator &cell,
unsigned int &fn,/*face number*/
FullMatrix<double> &cell_matrix,
unsigned int &i_dir,
unsigned int &g);

virtual void integrate_reflective_boundary_linear_form
(const std_cxx11::shared_ptr<FEFaceValues<dim> > fvf,
typename DoFHandler<dim>::active_cell_iterator &cell,
unsigned int &fn,/*face number*/
std::vector<Vector<double> > &cell_rhses,
unsigned int &i_dir,
unsigned int &g);

virtual void integrate_interface_bilinear_form
(const std_cxx11::shared_ptr<FEFaceValues<dim> > fvf,
const std_cxx11::shared_ptr<FEFaceValues<dim> > fvf_nei,
typename DoFHandler<dim>::active_cell_iterator &cell,
typename DoFHandler<dim>::cell_iterator &neigh,/*cell iterator for cell*/
unsigned int &fn,/*concerning face number in local cell*/
unsigned int &i_dir,
unsigned int &g,
FullMatrix<double> &vp_up,
FullMatrix<double> &vp_un,
FullMatrix<double> &vn_up,
FullMatrix<double> &vn_un);

virtual void generate_moments ();
virtual void postprocess ();
virtual void generate_ho_rhs ();
virtual void generate_ho_fixed_source ();

private:
void setup_system ();
void generate_globally_refined_grid ();
Expand Down
35 changes: 25 additions & 10 deletions src/common/bart_tools.cc
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,25 @@
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/mpi.h>

#include "bart_builder.h"
#include "../transport/even_parity.h"
#include "../equations/even_parity.h"
#include "../aqdata/aq_lsgc"

template <int dim>
FE_Poly<TensorProductPolynomials<dim>,dim,dim>* build_finite_element
(ParameterHandler &prm)
{
std::string discretization = prm.get ("spatial discretization");
AssertThrow (discretization!="none",
ExcMessage("valid spatial discretizations are dfem and cfem"))
if (discretization=="dfem")
return new FE_DGQ<dim> (p_order);
else
return new FE_Q<dim> (p_order);
}

template <int dim>
std_cxx11::shared_ptr<MeshGenerator<dim> > build_mesh (ParameterHandler &prm)
{
Expand All @@ -22,10 +37,10 @@ std_cxx11::shared_ptr<MaterialProperties> build_material (ParameterHandler &prm)
}

template <int dim>
std_cxx11::shared_ptr<IterationBase<dim> >
build_transport_iteration (ParameterHandler &prm,
const std_cxx11::shared_ptr<MeshGenerator<dim> > msh_ptr,
const std_cxx11::shared_ptr<AQBase<dim> > aqd_ptr)
std_cxx11::shared_ptr<IterationBase<dim> > build_transport_iteration
(ParameterHandler &prm,
const std_cxx11::shared_ptr<MeshGenerator<dim> > msh_ptr,
const std_cxx11::shared_ptr<AQBase<dim> > aqd_ptr)
{
// in future development this builder will be like other two
std_cxx11::shared_ptr<IterationBase<dim> > iteration_class =
Expand All @@ -36,11 +51,11 @@ build_transport_iteration (ParameterHandler &prm,
}

template <int dim>
std_cxx11::shared_ptr<TransportBase<dim> >
build_transport_model (ParameterHandler &prm,
const std_cxx11::shared_ptr<MeshGenerator<dim> > msh_ptr,
const std_cxx11::shared_ptr<AQBase<dim> > aqd_ptr,
const std_cxx11::shared_ptr<MaterialProperties> mat_ptr)
std_cxx11::shared_ptr<TransportBase<dim> > build_transport_model
(ParameterHandler &prm,
const std_cxx11::shared_ptr<MeshGenerator<dim> > msh_ptr,
const std_cxx11::shared_ptr<AQBase<dim> > aqd_ptr,
const std_cxx11::shared_ptr<MaterialProperties> mat_ptr)
{
std::string transport_model_name = prm.get ("transport model");
AssertThrow (transport_model_name!="none",
Expand Down
7 changes: 6 additions & 1 deletion src/common/bart_tools.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,14 +2,19 @@
#define __bart_builder_h__

#include <deal.II/base/parameter_handler.h>
#include <deal.II/fe/fe_poly.h>

#include <string>
#include <map>

#include "../transport/transport_base.h"
#include "../equations/equation_base.h"

using namespace dealii;

template <int dim>
FE_Poly<TensorProductPolynomials<dim>,dim,dim>*
build_finite_element (ParameterHandler &prm);

template <int dim>
std_cxx11::shared_ptr<MeshGenerator<dim> >
build_mesh (ParameterHandler &prm);
Expand Down
56 changes: 18 additions & 38 deletions src/iteration/iteration_base.cc
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
#include <deal.II/base/index_set.h>
#include <deal.II/fe/fe_values.h>

#include <boost/algorithm/string.hpp>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/grid/cell_id.h>

Expand All @@ -9,7 +8,7 @@

#include <algorithm>

#include "transport_base.h"
#include "equation_base.h"
#include "../aqdata/aq_base.h"
#include "../aqdata/aq_lsgc.h"

Expand All @@ -19,7 +18,8 @@ template <int dim>
IterationBase<dim>::IterationBase
(ParameterHandler &prm,
const std_cxx11::shared_ptr<MeshGenerator<dim> > msh_ptr,
const std_cxx11::shared_ptr<AQBase<dim> > aqd_ptr)
const std_cxx11::shared_ptr<AQBase<dim> > aqd_ptr,
const IndexSet local_dofs)
:
err_k_tol(1.0e-6),
err_phi_tol(1.0e-7),
Expand All @@ -39,6 +39,7 @@ pcout(std::cout,
is_cell_at_bd,
is_cell_at_ref_bd);
trm_ptr = build_transport_model ();
initialize_intermediate_vectors (local_dofs);
sflx_proc.resize (n_group);
sflx_proc_prev_gen.resize (n_group);
}
Expand Down Expand Up @@ -124,8 +125,9 @@ void IterationBase<dim>::power_iteration
{
ct += 1;
update_ho_moments_in_fiss ();
scale_fiss_transfer_matrices ();
generate_ho_fixed_source ();
trm_ptr->scale_fiss_transfer_matrices ();//???
trm_ptr->generate_ho_fixed_source (vec_ho_fixed_rhs,
sflx_proc_prev_gen);
source_iteration ();
update_fiss_source_keff ();
err_phi = estimate_phi_diff (vec_ho_sflx, vec_ho_sflx_prev_gen);
Expand All @@ -138,19 +140,24 @@ void IterationBase<dim>::power_iteration
}

template <int dim>
void IterationBase<dim>::source_iteration ()
void IterationBase<dim>::source_iteration
(std::vector<PETScWrappers::MPI::SparseMatrix*> &vec_ho_sys,
std::vector<PETScWrappers::MPI::Vector*> &vec_aflx,
std::vector<PETScWrappers::MPI::Vector*> &vec_ho_rhs,
std::vector<PETScWrappers::MPI::Vector*> &vec_ho_fixed_rhs,
)
{
unsigned int ct = 0;
double err_phi = 1.0;
double err_phi_old;
while (err_phi>err_phi_tol)
{
ct += 1;
generate_ho_rhs ();
trm_ptr->generate_ho_rhs ();
sol_ptr->ho_solve (vec_ho_sys,
vec_aflx,
vec_ho_rhs);
generate_moments ();
trm_ptr->generate_moments ();
err_phi_old = err_phi;
err_phi = estimate_phi_diff (vec_ho_sflx, vec_ho_sflx_old);
double spectral_radius = err_phi / err_phi_old;
Expand All @@ -166,33 +173,6 @@ void IterationBase<dim>::postprocess ()
{// do nothing in the base class
}

template <int dim>
double IterationBase<dim>::estimate_fiss_source (std::vector<Vector<double> > &phis_this_process)
{
double fiss_source = 0.0;
for (unsigned int ic=0; ic<local_cells.size(); ++ic)
{
typename DoFHandler<dim>::active_cell_iterator cell = local_cells[ic];
std::vector<std::vector<double> > local_phis (n_group,
std::vector<double> (n_q));
unsigned int material_id = cell->material_id ();
if (is_material_fissile[material_id])
{
fv->reinit (cell);
for (unsigned int g=0; g<n_group; ++g)
fv->get_function_values (phis_this_process[g],
local_phis[g]);
for (unsigned int qi=0; qi<n_q; ++qi)
for (unsigned int g=0; g<n_group; ++g)
fiss_source += (all_nusigf[material_id][g] *
local_phis[g][qi] *
fv->JxW(qi));
}
}
double global_fiss_source = Utilities::MPI::sum (fiss_source, mpi_communicator);
return global_fiss_source;
}

template <int dim>
double IterationBase<dim>::estimate_k (double &fiss_source,
double &fiss_source_prev_gen,
Expand All @@ -211,8 +191,8 @@ double IterationBase<dim>::estimate_phi_diff
double err = 0.0;
for (unsigned int i=0; i<phis_newer.size (); ++i)
{
PETScWrappers::MPI::Vector dif = *(phis_newer)[i];
dif -= *(phis_older)[i];
PETScWrappers::MPI::Vector dif = *phis_newer[i];
dif -= *phis_older[i];
err = std::max (err, dif.l1_norm () / phis_newer[i]->l1_norm ());
}
return err;
Expand Down
Loading

0 comments on commit 1bdf3d8

Please sign in to comment.