diff --git a/src/iteration/eigen_base.cc b/src/iteration/eigen_base.cc index c273b939e..2a2ff403c 100644 --- a/src/iteration/eigen_base.cc +++ b/src/iteration/eigen_base.cc @@ -15,10 +15,15 @@ EigenBase::~EigenBase () } template -EigenBase::do_iterations () +EigenBase::do_iterations +(ParameterHandler &prm, + std_cxx11::shared_ptr > msh_ptr, + std_cxx11::shared_ptr > aqd_ptr, + std_cxx11::shared_ptr mat_ptr, + std::vector sys_mats) { - this->initialize_equations (); - eigen_iterations (); + this->initialize_equations (prm, msh_ptr, aqd_ptr, mat_ptr); + eigen_iterations (msh_ptr, aqd_ptr, mat_ptr, sys_mats); } template diff --git a/src/iteration/in_group_base.h b/src/iteration/in_group_base.h index e62345c23..dbe9971df 100644 --- a/src/iteration/in_group_base.h +++ b/src/iteration/in_group_base.h @@ -1,6 +1,7 @@ #ifndef __in_group_base_h__ #define __in_group_base_h__ +template class InGroupBase : public IterationBase { public: diff --git a/src/iteration/iteration_base.cc b/src/iteration/iteration_base.cc index 8171e3438..56549720c 100644 --- a/src/iteration/iteration_base.cc +++ b/src/iteration/iteration_base.cc @@ -72,8 +72,15 @@ double IterationBase::estimate_phi_diff } template -void IterationBase::initialize_equations () +void IterationBase::initialize_equations +(ParameterHandler &prm, + std_cxx11::shared_ptr > msh_ptr, + std_cxx11::shared_ptr > aqd_ptr, + std_cxx11::shared_ptr mat_ptr) { + tra_ptr = build_transport_model (prm, msh_ptr, aqd_ptr, mat_ptr); + if (do_nda) + nda_ptr = build_nda (prm, msh_ptr, aqd_ptr, mat_ptr); } template class IterationBase<2>; diff --git a/src/iteration/iteration_base.h b/src/iteration/iteration_base.h index 2788cb9d7..14d3d7735 100644 --- a/src/iteration/iteration_base.h +++ b/src/iteration/iteration_base.h @@ -59,6 +59,9 @@ class IterationBase void initialize_equations (); + std_cxx11::shared_ptr > tra_ptr; + std_cxx11::shared_ptr > nda_ptr; + double total_calculation_time; /**< total time for calculations including assembly of rhs*/ unsigned int ct_ho_iters; /**< HO iteration counts*/ unsigned int ct_nda_iters; /**< NDA iteration counts*/ diff --git a/src/iteration/mg_base.cc b/src/iteration/mg_base.cc index e0cb793a8..53522f5b4 100644 --- a/src/iteration/mg_base.cc +++ b/src/iteration/mg_base.cc @@ -15,20 +15,46 @@ void MGBase::do_iterations (std::vector &sys_mats, std::vector &sys_rhses) { - this->initialize_equations (); - multigroup_iterations (); + this->initialize_equations (prm, msh_ptr, aqd_ptr, mat_ptr); + mg_iterations (msh_ptr, aqd_ptr, mat_ptr); } template -void MGBase::multigroup_iterations +void MGBase::mg_iterations (std::vector &sys_mats, std::vector &sys_rhses) {// this function needs to be overridden if JFNK is desired + + /* for (unsigned int g=0; gsolve_in_group (sys_mats, g) } + */ + // GS + /* + for (unsigned int g=0; gsolve_in_group (sys_mats,vec_aflx,sys_rhses) + } + */ + // Jacobi + /* + for (unsigned int g=0; gsolve_in_group (sys_mats,vec_aflx,sys_rhses) + */ +} + +template +void MGBase::generate_system_matrices +(std::vector &sys_mats, + std_cxx11::shared_ptr > equ_ptr) +{ + equ_ptr->assemble_system (); } template diff --git a/src/iteration/mg_base.h b/src/iteration/mg_base.h index 8d2a8656a..8b615de6f 100644 --- a/src/iteration/mg_base.h +++ b/src/iteration/mg_base.h @@ -10,6 +10,7 @@ class MGBase : public IterationBase MGBase (); virtual ~MGBase (); + virtual void mg_iterations (); virtual void generate_system_matrices (std::vector &sys_mats); virtual void generate_group_rhses diff --git a/src/source_iteration.cc b/src/source_iteration.cc new file mode 100644 index 000000000..e69de29bb diff --git a/src/source_iteration.h b/src/source_iteration.h new file mode 100644 index 000000000..e69de29bb diff --git a/src/transport/equation_base.cc b/src/transport/equation_base.cc index 55548a965..824359a43 100644 --- a/src/transport/equation_base.cc +++ b/src/transport/equation_base.cc @@ -34,7 +34,6 @@ p_order(prm.get_integer("finite element polynomial degree")), nda_quadrature_order(p_order+3) //this is hard coded { this->process_input (msh_ptr, aqd_ptr, mat_ptr); - } template @@ -57,7 +56,9 @@ void EquationBase::process_input // aq data related { - n_total_ho_vars = aqd_ptr->get_n_total_ho_vars (); + // note that n_total_vars will be have to be re-init + // in derived class of if it's for NDA + n_total_vars = aqd_ptr->get_n_total_ho_vars (); n_azi = aqd_ptr->get_sn_order (); n_dir = aqd_ptr->get_n_dir (); component_index = aqd_ptr->get_component_index_map (); @@ -137,6 +138,8 @@ void EquationBase::initialize_assembly_related_objects n_q = q_rule->size(); n_qf = qf_rule->size(); + // the following section will be for NDA soly + /* if (do_nda) { qc_rule = std_cxx11::shared_ptr > (new QGauss (nda_quadrature_order)); @@ -155,6 +158,7 @@ void EquationBase::initialize_assembly_related_objects n_qc = qc_rule->size (); n_qfc = qfc_rule->size (); } + */ local_dof_indices.resize (dofs_per_cell); neigh_dof_indices.resize (dofs_per_cell); @@ -162,7 +166,8 @@ void EquationBase::initialize_assembly_related_objects template void EquationBase::assemble_volume_boundary -(std::vector::active_cell_iterator> &local_cells, +(std::vector &sys_mats, + std::vector::active_cell_iterator> &local_cells, std::vector &is_cell_at_bd) { // volumetric pre-assembly matrices @@ -172,9 +177,6 @@ void EquationBase::assemble_volume_boundary std::vector > collision_at_qp (n_q, FullMatrix(dofs_per_cell, dofs_per_cell)); - for (unsigned int i=0; i (n_q, dofs_per_cell)); - // this sector is for pre-assembling streaming and collision matrices at quadrature // points { @@ -183,7 +185,7 @@ void EquationBase::assemble_volume_boundary pre_assemble_cell_matrices (fv, cell, streaming_at_qp, collision_at_qp); } - for (unsigned int k=0; k::assemble_volume_boundary local_mat, g, i_dir); } - - if (k==0) - for (unsigned int qi=0; qishape_value (i,qi) * fv->JxW (qi); - - vec_ho_sys[k]->add (local_dof_indices, + sys_mats[k]->add (local_dof_indices, local_dof_indices, local_mat); } - vec_ho_sys[k]->compress (VectorOperation::add); + sys_mats[k]->compress (VectorOperation::add); }// components } @@ -295,14 +291,15 @@ void EquationBase::integrate_reflective_boundary_linear_form */ template void EquationBase::assemble_interface -(std::vector::active_cell_iterator> &local_cells) +(std::vector &sys_mats, + std::vector::active_cell_iterator> &local_cells) { FullMatrix vi_ui (dofs_per_cell, dofs_per_cell); FullMatrix vi_ue (dofs_per_cell, dofs_per_cell); FullMatrix ve_ui (dofs_per_cell, dofs_per_cell); FullMatrix ve_ue (dofs_per_cell, dofs_per_cell); - for (unsigned int k=0; k::assemble_interface fn, vi_ui, vi_ue, ve_ui, ve_ue, g, i_dir/*specific component*/); - vec_ho_sys[k]->add (local_dof_indices, - local_dof_indices, - vi_ui); - - vec_ho_sys[k]->add (local_dof_indices, - neigh_dof_indices, - vi_ue); - - vec_ho_sys[k]->add (neigh_dof_indices, - local_dof_indices, - ve_ui); - - vec_ho_sys[k]->add (neigh_dof_indices, - neigh_dof_indices, - ve_ue); + sys_mats[k]->add (local_dof_indices, + local_dof_indices, + vi_ui); + + sys_mats[k]->add (local_dof_indices, + neigh_dof_indices, + vi_ue); + + sys_mats[k]->add (neigh_dof_indices, + local_dof_indices, + ve_ui); + + sys_mats[k]->add (neigh_dof_indices, + neigh_dof_indices, + ve_ue); }// target faces } - vec_ho_sys[k]->compress(VectorOperation::add); + sys_mats[k]->compress(VectorOperation::add); }// component } @@ -380,22 +377,21 @@ void EquationBase::integrate_interface_bilinear_form template void EquationBase::generate_moments -(std::vector &vec_ho_sflx, - std::vector &vec_ho_sflx_old, - std::vector > &sflx_proc) +(std::vector &vec_aflx, + std::vector > &sflx_proc, + std::vector > &sflx_proc_old) { // PETSc type vectors live in BartDriver // TODO: only scalar flux is generated for now, future will be moments - AssertThrow(do_nda==false, - ExcMessage("Moments are generated only without NDA")); + //AssertThrow(do_nda==false, + // ExcMessage("Moments are generated only without NDA")); if (!do_nda) for (unsigned int g=0; gadd (wi[i_dir], *vec_aflx[get_component_index(i_dir, g)]); - sflx_proc[g] = *vec_ho_sflx[g]; + sflx_proc[g].add (wi[i_dir], *vec_aflx[get_component_index(i_dir, g)]); } } @@ -420,20 +416,21 @@ void EquationBase::scale_fiss_transfer_matrices (double keff) } template -void EquationBase::generate_ho_rhs -(std::vector &vec_ho_rhs, - std::vector &vec_ho_fixed_rhs, - std::vector &sflx_this_proc) +void EquationBase::generate_rhs +(std::vector &vec_fixed_rhs, + std::vector > &sflx_proc, + std::vector > &sflx_proc_old) { } -void EquationBase::generate_ho_fixed_source +void EquationBase::generate_fixed_source (std::vector &vec_ho_fixed_rhs, std::vector > &sflx_this_proc) template double EquationBase::estimate_fiss_source -(std::vector > &phis_this_process) +(std::vector > &phis_this_process, + std::vector::active_cell_iterator> &local_cells) { // first, estimate local fission source double fiss_source = 0.0; diff --git a/src/transport/equation_base.h b/src/transport/equation_base.h index afbb1030e..27d520e10 100644 --- a/src/transport/equation_base.h +++ b/src/transport/equation_base.h @@ -45,10 +45,17 @@ template class EquationBase { public: + // For SN and NDA (NDA still needs quadrature to calculate corrections) EquationBase (ParameterHandler &prm, const std_cxx11::shared_ptr > msh_ptr, const std_cxx11::shared_ptr > aqd_ptr, const std_cxx11::shared_ptr mat_ptr) + + // For diffusion (future work for someone) + EquationBase (ParameterHandler &prm, + const std_cxx11::shared_ptr > msh_ptr, + const std_cxx11::shared_ptr mat_ptr) + virtual ~EquationBase (); void run (); @@ -69,7 +76,8 @@ class EquationBase std::vector > &collision_at_qp); virtual void integrate_cell_bilinear_form - (const std_cxx11::shared_ptr > fv, + (std::vector &sys_mats, + const std_cxx11::shared_ptr > fv, typename DoFHandler::active_cell_iterator &cell, FullMatrix &cell_matrix, std::vector > > &streaming_at_qp, @@ -78,7 +86,8 @@ class EquationBase const unsigned int &i_dir=0); virtual void integrate_boundary_bilinear_form - (const std_cxx11::shared_ptr > fvf, + (std::vector &sys_mats, + const std_cxx11::shared_ptr > fvf, typename DoFHandler::active_cell_iterator &cell, unsigned int &fn,/*face number*/ FullMatrix &cell_matrix, @@ -94,7 +103,8 @@ class EquationBase const unsigned int &i_dir=0); virtual void integrate_interface_bilinear_form - (const std_cxx11::shared_ptr > fvf, + (std::vector &sys_mats, + const std_cxx11::shared_ptr > fvf, const std_cxx11::shared_ptr > fvf_nei, typename DoFHandler::active_cell_iterator &cell, typename DoFHandler::cell_iterator &neigh,/*cell iterator for cell*/ @@ -107,15 +117,20 @@ class EquationBase FullMatrix &vn_un); virtual void generate_moments - (std::vector &vec_ho_sflx, - std::vector &vec_ho_sflx_old, - std::vector > &sflx_proc); + (std::vector*> &vec_ho_sflx, + std::vector*> &vec_ho_sflx_old, + std::vector*> &sflx_proc); virtual void postprocess (); virtual void generate_ho_rhs (); virtual void generate_ho_fixed_source (std::vector &vec_ho_fixed_rhs, std::vector > &sflx_this_proc); + // override this three functions in derived classes + virtual unsigned int get_component_index (unsigned int incident_angle_index, unsigned int g); + virtual unsigned int get_component_direction (unsigned int comp_ind); + virtual unsigned int get_component_group (unsigned int comp_ind); + private: void setup_system (); void generate_globally_refined_grid (); @@ -156,10 +171,6 @@ class EquationBase std::string aq_name; protected: - unsigned int get_component_index (unsigned int incident_angle_index, unsigned int g); - unsigned int get_component_direction (unsigned int comp_ind); - unsigned int get_component_group (unsigned int comp_ind); - unsigned int get_reflective_direction_index (unsigned int boundary_id, unsigned int incident_angle_index); @@ -192,7 +203,7 @@ class EquationBase unsigned int n_dir; unsigned int n_azi; - unsigned int n_total_ho_vars; + unsigned int n_total_vars; unsigned int n_group; unsigned int n_material; unsigned int p_order; diff --git a/src/transport/even_parity.h b/src/transport/even_parity.h index 4ac67bf68..14672ab7e 100644 --- a/src/transport/even_parity.h +++ b/src/transport/even_parity.h @@ -49,11 +49,11 @@ class EvenParity : public EquationBase const unsigned int &g, const unsigned int &i_dir); - void generate_ho_fixed_source + void generate_fixed_source (std::vector &vec_ho_fixed_rhs, std::vector > &sflx_this_proc); - void generate_ho_rhs + void generate_rhs (std::vector &vec_ho_rhs, std::vector &vec_ho_fixed_rhs, std::vector &sflx_this_proc);