All the changes made to make SlaybaughLab#10 work (in progress)

src/iteration/eigen_base.cc

@@ -1,11 +1,13 @@
 #include "eigen_base.h"
+#include "mg_base.h"
 
 template <int dim>
-EigenBase<dim>::EigenBase () : IterationBase<dim> (),
+EigenBase<dim>::EigenBase (ParameterHandler &prm) : IterationBase<dim> (prm),
 err_k_tol(1.0e-6),
 err_phi_tol(1.0e-5),
 keff(1.0)
 {
+  mg_ptr = build_mg_iterations (prm);
 }
 
 template <int dim>
@@ -15,39 +17,48 @@ EigenBase<dim>::~EigenBase ()
 
 template <int dim>
 EigenBase<dim>::do_iterations
-(ParameterHandler &prm,
- std_cxx11::shared_ptr<MeshGenerator<dim> > msh_ptr,
- std_cxx11::shared_ptr<AQBase<dim> > aqd_ptr,
- std_cxx11::shared_ptr<MaterialProperties> mat_ptr,
- std::vector<PETScWrappers::MPI::SparseMatrix*> sys_mats)
+void MGBase<dim>::do_iterations
+(std::vector<Vector<double> > &sflx_proc,
+ std::vector<std_cxx11::shared_ptr<EquationBase<dim> > > &equ_ptrs)
 {
-  this->initialize_equations (prm, msh_ptr, aqd_ptr, mat_ptr);
-  eigen_iterations (msh_ptr, aqd_ptr, mat_ptr, sys_mats);
+  // assemble system matrices
+  for (unsigned int i=0; i<equ_ptrs.size(); ++i)
+    equ_ptrs[i]->assemble_bilinear_form ();
+
+  // initialize fission process
+  initialize_fiss_process (sflx_proc, equ_ptrs);
+
+  // perform eigenvalue iterations
+  eigen_iterations (sflx_proc, equ_ptrs);
 }
 
 template <int dim>
-void EigenBase<dim>::generate_system_matrices
-(std::vector<PETScWrappers::MPI::SparseMatrix*> sys_mats)
+void EigenBase<dim>::initialize_fiss_process
+(std::vector<Vector<double> > &sflx_proc,
+ std::vector<std_cxx11::shared_ptr<EquationBase<dim> > > &equ_ptrs)
 {
-  // EquationBase will be instantiated in InGroupBase
-  mgs_ptr->generate_system_matrices (sys_mats);
+  for (unsigned int g=0; g<n_group; ++g)
+    sflx_proc[g] = 1.0;
+
+  fission_source = equ_ptrs[0]->estimate_fiss_source (this->sflx_proc);
 }
 
+// Override this function to do specific eigenvalue iteration as desired
 template <int dim>
-void EigenBase<dim>::initialize_fiss_process ()
+void EigenBase<dim>::eigen_iterations
+(std::vector<Vector<double> > &sflx_proc,
+ std::vector<std_cxx11::shared_ptr<EquationBase<dim> > > &equ_ptrs)
 {
-  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>
-void EigenBase<dim>::update_fiss_source_keff ()
+void EigenBase<dim>::update_fiss_source_keff
+(std::vector<Vector<double> > &sflx_proc,
+ std::vector<std_cxx11::shared_ptr<EquationBase<dim> > > &equ_ptrs)
 {
-  keff_prev_gen = keff;
-  fission_source_prev_gen = fission_source;
-  fission_source = trm_ptr->estimate_fiss_source (sflx_proc);
+  keff_prev = keff;
+  fission_source_prev = fission_source;
+  fission_source = equ_ptrs[0]->estimate_fiss_source (sflx_proc);
   keff = estimate_k (fission_source, fission_source_prev_gen, keff_prev_gen);
 }
 
@@ -60,17 +71,11 @@ double EigenBase<dim>::estimate_k (double &fiss_source,
 }
 
 template <int dim>
-double EigenBase<dim>::estimate_k_err (double &k, double &k_prev)
+double EigenBase<dim>::estimate_k_diff (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 <int dim>
 double EigenBase<dim>::get_keff ()
 {

src/iteration/eigen_base.h

@@ -1,14 +1,17 @@
 #ifndef __eigen_base_h__
 #define __eigen_base_h__
 
+#include "iteration_base.h"
+
 template <int dim>
 class EigenBase : public IterationBase
 {
 public:
-  EigenBase ();
+  EigenBase (ParameterHandler &prm);
   virtual ~EigenBase ();
 
   virtual void do_iterations ();
+
   virtual void eigen_iterations ();
 
   double get_keff ();
@@ -21,6 +24,8 @@ class EigenBase : public IterationBase
   const double err_k_tol;
   const double err_phi_tol;
 
+  std_cxx11::shared_ptr<MGBase<dim> > mg_ptr;
+
   double keff;
   double keff_prev;
 

src/iteration/in_group_base.cc

@@ -1,12 +1,12 @@
 #include "in_group_base.h"
-#include "../common/preconditioner_solver.h"
 
 template <int dim>
 InGroupBase<dim>::InGroupBase ()
 :
 IterationBase<dim> (),
 err_phi_tol(1.0e-6)
 {
+  sflx_proc_old.resize (1);
 }
 
 template <int dim>
@@ -16,7 +16,8 @@ InGroupBase<dim>::~InGroupBase ()
 
 template <int dim>
 InGroupBase<dim>::solve_in_group
-(std_cxx11::shared_ptr<EquationBase<dim> > equ_ptr,
+(std::vector<Vector<double> > &sflx_proc,
+ std_cxx11::shared_ptr<EquationBase<dim> > equ_ptrs,
  unsigned int &g)
 {
 }
@@ -30,11 +31,22 @@ InGroupBase<dim> ()
 
 template <int dim>
 SourceIteration<dim>::solve_in_group
-(std_cxx11::shared_ptr<EquationBase<dim> > equ_ptr,
+(std::vector<Vector<double> > &sflx_proc,
+ std_cxx11::shared_ptr<EquationBase<dim> > equ_ptrs,
  unsigned int &g)
 {
   double err = 1.0;
-  while (err<this->err_phi_tol)
+  while (err>this->err_phi_tol)
+  {
+    // generate rhs for group g
+    equ_ptrs[0]->assemble_linear_form (sflx_proc, g);
+    // solve all the directions in group g
+    equ_ptrs[0]->solve_in_group (g);
+    // generate moments
+    equ_ptrs[0]->generate (sflx_proc[g], sflx_proc_old, g);
+    // calculate the difference of moments for convergence check
+    err = this->estimate_phi_diff (sflx_proc[g], sflx_proc_old);
+  }
 }
 
 template class InGroupBase<2>;

src/iteration/in_group_base.h

@@ -1,35 +1,37 @@
 #ifndef __in_group_base_h__
 #define __in_group_base_h__
 
+#include "in_group_base.h"
+
 template <int dim>
 class InGroupBase : public IterationBase<dim>
 {
 public:
-  InGroupBase ();
+  InGroupBase (ParameterHandler &prm);
   virtual ~InGroupBase();
 
   // has to be provided
   virtual void solve_in_group
-  (std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells,
-   std::vector<bool> &is_cell_at
+  (std::vector<Vector<double> > &sflx_proc,
    std_cxx11::shared_ptr<EquationBase<dim> > equ_ptrs,
    unsigned int &g);
 
 protected:
   const double err_phi_tol;
 
-  std::vector<Vector<double> > sflx_proc_old;
+  Vector<double> sflx_proc_old;
 };
 
 template <int dim>
 class SourceIteration : public InGroupBase<dim>
 {
 public:
-  SourceIteration ();
+  SourceIteration (ParameterHandler &prm);
   ~SourceIteration ();
 
   void solve_in_group
-  (std_cxx11::shared_ptr<EquationBase<dim> > equ_ptrs,
+  (std::vector<Vector<double> > &sflx_proc,
+   std_cxx11::shared_ptr<EquationBase<dim> > equ_ptrs,
    unsigned int &g);
 };
 

src/iteration/iteration_base.cc

@@ -1,7 +1,7 @@
 #include "iteration_base.h"
 
 template <int dim>
-IterationBase<dim>::IterationBase ()
+IterationBase<dim>::IterationBase (ParameterHandler &prm)
 {
 }
 
@@ -10,6 +10,13 @@ IterationBase<dim>::~IterationBase ()
 {
 }
 
+template <int dim>
+void IterationBase<dim>::do_iterations
+(std::vector<Vector<double> > &sflx_proc,
+ std::vector<std_cxx11::shared_ptr<EquationBase<dim> > > &equ_ptrs)
+{
+}
+
 template <int dim>
 double IterationBase<dim>::estimate_phi_diff
 (std::vector<PETScWrappers::MPI::Vector*> &phis_newer,

src/iteration/iteration_base.h

@@ -1,18 +1,20 @@
 #ifndef __iteration_base_h__
 #define __iteration_base_h__
 
+#include <deal.II/base/parameter_handler.h>
+
 using namespace dealii;
 
 template <int dim>
 class IterationBase
 {
 public:
-  IterationBase ();
+  IterationBase (ParameterHandler &prm);
   virtual ~IterationBase ();
 
-  virtual void do_iterations ();
-  virtual void assemble_system
-  (std::vector<PETScWrappers::MPI::SparseMatrix*> sys_mats);
+  virtual void do_iterations
+  (std::vector<Vector<double> > &sflx_proc,
+   std::vector<std_cxx11::shared_ptr<EquationBase<dim> > > &equ_ptrs);
 
 protected:
   /** \brief Function to measure the relative difference between two sets of PETSc

src/iteration/iterations.cc

@@ -43,6 +43,8 @@ template <int dim>
 void Iterations<dim>::initialize_system_matrices_vectors
 (SparsityPatternType &dsp, IndexSet &local_dofs)
 {
+  // each equation pointer contains system matrices and vectors from PETSc, per se
+  // This function is to tell the shapes and sizes of those vectors and matrices
   for (unsigned int i=0; i<equ_ptrs.size(); ++i)
     equ_ptrs[i]->initialize_system_matrices_vectors (dsp, local_dofs);
 }
@@ -53,25 +55,22 @@ void Iterations<dim>::solve_problems (std::vector<Vector<double> > &sflx_proc)
   if (is_eigen_problem)
   {
     std_cxx11::shared_ptr<EigenBase<dim> > pro_ptr = build_eigen_problem (prm);
-    pro_ptr->do_iterations (local_cells, is_cell_at_bd, sflx_proc, equ_ptrs);
+    pro_ptr->do_iterations (sflx_proc, equ_ptrs);
     pro_ptr->get_keff (keff);
   }
   else
   {
     std_cxx11::shared_ptr<MGBase<dim> > pro_ptr = build_mg_problem (prm);
-    pro_ptr->do_iterations (local_cells, is_cell_at_bd, sflx_proc, equ_ptrs);
+    pro_ptr->do_iterations (sflx_proc, equ_ptrs);
   }
 }
 
-template <int dim>
-void Iterations<dim>
-
 template <int dim>
 void Iterations<dim>::get_keff (double &k)
 {
   AssertThrow (is_eigen_problem,
                ExcMessage("Only eigen problems have keff"));
-  k = this->keff;
+  k = keff;
 }
 
 // explicit instantiation to avoid linking error

src/iteration/mg_base.cc

@@ -1,10 +1,11 @@
 #include "mg_base.h"
 
 template <int dim>
-MGBase<dim>::MGBase ()
-: IterationBase<dim> (),
+MGBase<dim>::MGBase (ParameterHandler &prm)
+: IterationBase<dim> (prm),
 err_phi_tol(1.0e-5)
 {
+  ig_ptr = build_ig_iteration (prm);
 }
 
 template <int dim>
@@ -14,29 +15,27 @@ MGBase<dim>::~MGBase ()
 
 template <int dim>
 void MGBase<dim>::do_iterations
-(std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells,
- std::vector<bool> &is_cell_at_bd,
+(std::vector<Vector<double> > &sflx_proc,
  std::vector<std_cxx11::shared_ptr<EquationBase<dim> > > &equ_ptrs)
 {
-  // assemble bilinear form of transport equation
-  equ_ptrs[0]->assemble_bilinear_forms (local_cells, is_cell_at_bd);
-  // assemble NDA bilinear form if do_nda
-  if (do_nda)
-    equ_ptrs[1]->assemble_bilinear_forms (local_cells, is_cell_at_bd);
+  // assemble bilinear forms of available equations
+  for (unsigned int i=0; i<equ_ptrs.size(); ++i)
+    equ_ptrs[i]->assemble_bilinear_forms ();
+
   // multigroup iterations
-  mg_iterations (local_cells, is_cell_at_bd, equ_ptrs);
+  mg_iterations (sflx_proc, equ_ptrs);
 }
 
+// virtual function for all multigroup iteration method. It has to be overriden
+// per derived class of MGBase. If it's fixed source problem, it will be called
+// internally in do_iterations. Otherwise, it will be called externally in EigenBase
+// instances.
 template <int dim>
 void MGBase<dim>::mg_iterations
-(std::vector<typename DoFHandler<dim>::active_cell_iterator> &local_cells,
- std::vector<bool> &is_cell_at_bd,
+(std::vector<Vector<double> > &sflx_proc,
  std::vector<std_cxx11::shared_ptr<EquationBase<dim> > > &equ_ptrs)
 {// this function needs to be overridden if JFNK is desired
   // by default, we give out Jacobi iteration scheme
-  for (unsigned int g=0; g<n_group; <#increment#>) {
-    <#statements#>
-  }
   /*
   for (unsigned int g=0; g<n_group; ++g)
   {