Skip to content

Commit

Permalink
Merge pull request #2 from cpraveen/master
Browse files Browse the repository at this point in the history 
Changed cell_number function. Removed MOOD and some other things which w...
  • Loading branch information
@juanpablogallego
juanpablogallego committed Oct 9, 2014
2 parents 8d2b56d + 24d3fb4 commit 2934889
Show file tree
Hide file tree
Showing 5 changed files with 8 additions and 829 deletions.
142 changes: 4 additions & 138 deletions src_mpi/claw.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,10 +79,6 @@ ConservationLaw<dim>::ConservationLaw (const char *input_filename,
verbose_cout (std::cout, false)
{
read_parameters (input_filename);

// For MOOD method compute degree reduction matrices
if(parameters.solver == Parameters::Solver::mood)
compute_reduction_matrices ();
}

//------------------------------------------------------------------------------
Expand All @@ -100,22 +96,6 @@ dh_cell (triangulation),
verbose_cout (std::cout, false)
{
read_parameters (input_filename);

// create map from dof index to total degree of basis function
index_to_degree.resize(fe.base_element(0).dofs_per_cell);
unsigned int c = 0;
if(dim==2)
{
for(unsigned int j=0; j<=degree; ++j)
for(unsigned int i=0; i<=degree-j; ++i)
{
index_to_degree[c++] = i+j;
}
}
else
{
AssertThrow(false, ExcMessage("Not implemented for dim=3"));
}
}

//------------------------------------------------------------------------------
Expand Down Expand Up @@ -305,21 +285,10 @@ void ConservationLaw<dim>::setup_system ()
system_matrix.reinit (sparsity_pattern);
}

// Allocate memory for MOOD variables
if(parameters.solver == Parameters::Solver::mood)
{
min_mood_var.reinit(triangulation.n_active_cells());
max_mood_var.reinit(triangulation.n_active_cells());
work1.reinit (dof_handler.n_dofs());
}

cell_degree.resize(triangulation.n_active_cells());
re_update.resize(triangulation.n_active_cells());
for(unsigned int i=0; i<triangulation.n_active_cells(); ++i)
{
cell_degree[i] = fe.degree;
re_update[i] = true;
}
unsigned int index=0;
for (typename DoFHandler<dim>::active_cell_iterator cell=dof_handler.begin_active();
cell!=dof_handler.end(); ++cell, ++index)
cell_number_map.insert (std::make_pair (std::pair<int,int>(cell->level(),cell->index()), index));

if(parameters.mapping_type == Parameters::AllParameters<dim>::cartesian)
compute_cartesian_mesh_size ();
Expand Down Expand Up @@ -573,7 +542,6 @@ ConservationLaw<dim>::compute_cell_average ()
for (; cell!=endc; ++cell)
{
unsigned int cell_no = cell_number(cell);
if(re_update[cell_no])
{
fe_values.reinit (cell);
fe_values.get_function_values (current_solution, solution_values);
Expand Down Expand Up @@ -760,105 +728,6 @@ void ConservationLaw<dim>::iterate_explicit (IntegratorExplicit<dim>& integrator
}
}

//------------------------------------------------------------------------------
// Perform SSPRK step using MOOD method
//------------------------------------------------------------------------------
template <int dim>
void ConservationLaw<dim>::iterate_mood (IntegratorExplicit<dim>& integrator,
Vector<double>& newton_update,
double& res_norm0, double& res_norm)
{
std::pair<unsigned int, double> convergence;

// Loop for RK stages
for(unsigned int rk=0; rk<n_rk; ++rk)
{
std::cout << "RK stage " << rk+1 << std::endl;
std::cout << "\t Iter n_reduce n_re_update n_reset\n";

// set time in boundary condition
// NOTE: We need to check if this is time accurate.
for (unsigned int boundary_id=0; boundary_id<Parameters::AllParameters<dim>::max_n_boundaries;
++boundary_id)
{
parameters.boundary_conditions[boundary_id].values.set_time(elapsed_time);
}

compute_min_max_mood_var();

// iterate forward euler until DMP is satisfied
bool terminate = false;
unsigned int mood_iter = 0;
predictor = current_solution;
shock_indicator = 0;

while(!terminate)
{
assemble_system (integrator);

res_norm = right_hand_side.l2_norm();
if(rk == 0) res_norm0 = res_norm;

convergence = solve (newton_update, res_norm);

if(mood_iter == 0)
{
// In first iteration all cells must be updated
current_solution += newton_update;
work1 = current_solution;
}
else
{
// Update cells with re_update == true
std::vector<unsigned int> dof_indices(fe.dofs_per_cell);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for(; cell != endc; ++cell)
{
unsigned int c = cell_number(cell);
if(re_update[c] == true)
{
cell->get_dof_indices( dof_indices );
for(unsigned int i=0; i<fe.dofs_per_cell; ++i)
{
current_solution(dof_indices[i]) += newton_update(dof_indices[i]);
unsigned int base_i = fe.system_to_component_index(i).second;
if(index_to_degree[base_i] > cell_degree[c])
current_solution(dof_indices[i]) = 0.0;
work1(dof_indices[i]) = current_solution(dof_indices[i]);
}
}
}
}

compute_cell_average ();
unsigned int n_reduce=0, n_re_update=0, n_reset=0;
terminate = apply_mood (n_reduce, n_re_update, n_reset);

++mood_iter;
std::printf("%12d %12d %12d %12d\n", mood_iter, n_reduce, n_re_update, n_reset);
}

// Forward euler has been accepted.
// Now update rk stage.
current_solution = work1;

// current_solution = ark*old_solution + (1-ark)*current_solution
current_solution.sadd (1.0-ark[rk], ark[rk], old_solution);

// Reset degree and update flags
for(unsigned int i=0; i<triangulation.n_active_cells(); ++i)
{
cell_degree[i] = fe.degree;
re_update[i] = true;
}
compute_cell_average ();
apply_limiter();
if(parameters.pos_lim) apply_positivity_limiter ();
}
}

//------------------------------------------------------------------------------
// Perform one step of implicit scheme
//------------------------------------------------------------------------------
Expand Down Expand Up @@ -1041,9 +910,6 @@ void ConservationLaw<dim>::run ()
}
else if(parameters.solver == Parameters::Solver::mood)
{
IntegratorExplicit<dim> integrator_explicit (dof_handler);
setup_mesh_worker (integrator_explicit);
iterate_mood(integrator_explicit, newton_update, res_norm0, res_norm);
}
else
{
Expand Down
39 changes: 4 additions & 35 deletions src_mpi/claw.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,30 +120,12 @@ class ConservationLaw
void compute_angular_momentum ();
void compute_cell_average ();
void apply_limiter ();
void apply_limiter_TVB_Qk_deprecated ();
void apply_limiter_TVB_Qk ();
void apply_limiter_TVB_Pk ();
void apply_positivity_limiter ();
void compute_shock_indicator ();
void compute_shock_indicator_u2 ();
void compute_shock_indicator_kxrcf ();

void compute_reduction_matrices();
void compute_min_max_mood_var();
bool apply_mood(unsigned int&, unsigned int&, unsigned int&);
void reduce_degree(const typename dealii::DoFHandler<dim>::cell_iterator&,
const unsigned int,
dealii::FEValues<dim>&);
void reduce_degree_Pk(const typename dealii::DoFHandler<dim>::cell_iterator&,
const unsigned int,
dealii::FEValues<dim>&);
void reduce_degree_Qk(const typename dealii::DoFHandler<dim>::cell_iterator&,
const unsigned int,
dealii::FEValues<dim>&);
void get_mood_second_derivatives(const typename dealii::DoFHandler<dim>::cell_iterator &cell,
std::vector<double>& D2);
bool test_u2(const typename dealii::DoFHandler<dim>::cell_iterator &cell);

void compute_mu_shock ();
void shock_cell_term (DoFInfo& dinfo, CellInfo& info);
void shock_boundary_term (DoFInfo& dinfo, CellInfo& info);
Expand Down Expand Up @@ -221,6 +203,7 @@ class ConservationLaw
dealii::Vector<double> mu_shock;
dealii::Vector<double> shock_indicator;
dealii::Vector<double> jump_indicator;
std::map<std::pair<int,int>,unsigned int> cell_number_map;

double global_dt;
double elapsed_time;
Expand Down Expand Up @@ -252,18 +235,9 @@ class ConservationLaw
// distributing the degrees of freedom.
dealii::SparseMatrix<double> system_matrix;
dealii::SparsityPattern sparsity_pattern;

// MOOD data
std::vector<unsigned int> cell_degree;
std::vector<bool> re_update;
std::vector<dealii::FullMatrix<double> > Rmatrix;
dealii::Vector<double> min_mood_var, max_mood_var;

std::vector< dealii::Vector<double> > inv_mass_matrix;

// For FE_DGP, maps dof index to its degree
std::vector<unsigned int> index_to_degree;

Parameters::AllParameters<dim> parameters;
dealii::ConditionalOStream verbose_cout;

Expand Down Expand Up @@ -328,21 +302,16 @@ class ConservationLaw
// Given a cell iterator, return the cell number
template <typename ITERATOR>
inline
unsigned int cell_number (const ITERATOR &cell) const
unsigned int cell_number (const ITERATOR &cell)
{
dealii::DoFAccessor<dim,dealii::DoFHandler<dim>,false>
dof_accessor (&triangulation,
cell->level(),
cell->index(),
&dh_cell);
return dof_accessor.dof_index(0);
return cell_number_map[std::make_pair(cell->level(),cell->index())];
}

// If cell is active, return cell average.
// If cell is not active, return area average of child cells.
inline
void get_cell_average(const typename dealii::DoFHandler<dim>::cell_iterator& cell,
dealii::Vector<double>& avg) const
dealii::Vector<double>& avg)
{
if(cell->active())
{
Expand Down
20 changes: 0 additions & 20 deletions src_mpi/indicator.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,32 +20,12 @@ void ConservationLaw<dim>::compute_shock_indicator ()
{
shock_indicator = 1e20;
}
else if(parameters.shock_indicator_type == Parameters::Limiter::u2)
{
compute_shock_indicator_u2();
}
else
{
compute_shock_indicator_kxrcf();
}
}

//-----------------------------------------------------------------------------
template <int dim>
void ConservationLaw<dim>::compute_shock_indicator_u2 ()
{
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();

for(; cell != endc; ++cell)
{
bool u2 = test_u2(cell);
unsigned int c = cell_number(cell);
shock_indicator[c] = (u2) ? 0.0 : 1.0e20;
}
}

//-----------------------------------------------------------------------------
template <int dim>
void ConservationLaw<dim>::compute_shock_indicator_kxrcf ()
Expand Down
Loading

0 comments on commit 2934889

Please sign in to comment.