Skip to content

Commit

Permalink
Merge pull request #333 from asartori86/fix_examples
Browse files Browse the repository at this point in the history 
Fix examples
  • Loading branch information
@asartori86
asartori86 authored Jul 14, 2016
2 parents 9558556 + b35bc60 commit c6d0651
Show file tree
Hide file tree
Showing 7 changed files with 171 additions and 130 deletions.
3 changes: 2 additions & 1 deletion examples/dynamic_stokes/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,8 @@ SET(_main source/main.cc)
# Normally you shouldn't need to change anything below.
############################################################
# Declare all source files the target consists of:
file(GLOB _files source/*cc)
file(GLOB _files source/*cc
include/*.h)

# Don't compile the main file into the library
IF(NOT "${_main}" STREQUAL "")
Expand Down
63 changes: 25 additions & 38 deletions examples/dynamic_stokes/include/stokes_ida.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,6 @@
#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/block_linear_operator.h>

#include <deal2lkit/sundials_interface.h>
#include <deal2lkit/ida_interface.h>
#include <deal2lkit/parameter_acceptor.h>
#include <deal2lkit/parsed_grid_generator.h>
Expand All @@ -60,69 +59,57 @@ using namespace deal2lkit;

typedef TrilinosWrappers::MPI::BlockVector VEC;
template<int dim>
class Stokes : public SundialsInterface<VEC>, public ParameterAcceptor
class Stokes : public ParameterAcceptor
{
public:

Stokes (const MPI_Comm &comm);
Stokes (const MPI_Comm comm);

virtual void declare_parameters (ParameterHandler &prm);

void run ();

/*********************************************************
* Public interface from SundialsInterface
*********************************************************/
virtual shared_ptr<VEC>
shared_ptr<VEC>
create_new_vector() const;

/** Returns the number of degrees of freedom. Pure virtual function. */
virtual unsigned int n_dofs() const;
/** Returns the number of degrees of freedom.*/
unsigned int n_dofs() const;

/** This function is called at the end of each iteration step for
* the ode solver. Once again, the conversion between pointers and
* other forms of vectors need to be done inside the inheriting
* class. */
virtual void output_step(const double t,
const VEC &solution,
const VEC &solution_dot,
const unsigned int step_number,
const double h);
void output_step(const double t,
const VEC &solution,
const VEC &solution_dot,
const unsigned int step_number);

/** This function will check the behaviour of the solution. If it
* is converged or if it is becoming unstable the time integrator
* will be stopped. If the convergence is not achived the
* calculation will be continued. If necessary, it can also reset
* the time stepper. */
virtual bool solver_should_restart(const double t,
const unsigned int step_number,
const double h,
VEC &solution,
VEC &solution_dot);
bool solver_should_restart(const double t,
VEC &solution,
VEC &solution_dot);

/** For dae problems, we need a
residual function. */
virtual int residual(const double t,
const VEC &src_yy,
const VEC &src_yp,
VEC &dst);
int residual(const double t,
const VEC &src_yy,
const VEC &src_yp,
VEC &dst);

/** Setup Jacobian system and preconditioner. */
virtual int setup_jacobian(const double t,
const VEC &src_yy,
const VEC &src_yp,
const VEC &residual,
const double alpha);
int setup_jacobian(const double t,
const VEC &src_yy,
const VEC &src_yp,
const double alpha);


/** Inverse of the Jacobian vector product. */
virtual int solve_jacobian_system(const double t,
const VEC &y,
const VEC &y_dot,
const VEC &residual,
const double alpha,
const VEC &src,
VEC &dst) const;
int solve_jacobian_system(const VEC &src,
VEC &dst) const;



Expand All @@ -132,7 +119,7 @@ class Stokes : public SundialsInterface<VEC>, public ParameterAcceptor
corresponding variable is a
differential component, zero
otherwise. */
virtual VEC &differential_components() const;
VEC &differential_components() const;

private:
void refine_mesh ();
Expand All @@ -150,7 +137,7 @@ class Stokes : public SundialsInterface<VEC>, public ParameterAcceptor

void update_constraints(const double &t);

const MPI_Comm &comm;
const MPI_Comm comm;

unsigned int initial_global_refinement;
unsigned int max_time_iterations;
Expand Down Expand Up @@ -205,7 +192,7 @@ class Stokes : public SundialsInterface<VEC>, public ParameterAcceptor

ParsedDataOut<dim, dim> data_out;

IDAInterface<VEC> dae;
IDAInterface<VEC> ida;

IndexSet global_partitioning;
std::vector<IndexSet> partitioning;
Expand Down
98 changes: 66 additions & 32 deletions examples/dynamic_stokes/source/stokes_ida.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,17 +40,17 @@
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/fe/mapping_q.h>

#include <sundials/sundials_types.h>
#include <nvector/nvector_parallel.h>
#include <sundials/sundials_math.h>

#include <deal.II/lac/solver_gmres.h>
template <int dim>
Stokes<dim>::Stokes (const MPI_Comm &communicator)
Stokes<dim>::Stokes (const MPI_Comm communicator)
:
SundialsInterface<VEC> (communicator),
comm(communicator),
comm(Utilities::MPI::duplicate_communicator(communicator)),
pcout (std::cout,
(Utilities::MPI::this_mpi_process(comm)
== 0)),
Expand All @@ -77,7 +77,7 @@ Stokes<dim>::Stokes (const MPI_Comm &communicator)
dirichlet_dot("Dirichlet dot", 3, "u,u,p", "0=u"),

data_out("Output Parameters", "vtu"),
dae(*this)
ida("IDA Solver Parameters", comm)
{}

template <int dim>
Expand Down Expand Up @@ -152,7 +152,7 @@ void Stokes<dim>::setup_dofs (const bool &first_run)
dof_handler->distribute_dofs (*fe);
DoFRenumbering::component_wise (*dof_handler, sub_blocks);

mapping = SP(new MappingQ<dim>(2));
mapping = SP(new MappingQ<dim,dim>(2));

dofs_per_block.clear();
dofs_per_block.resize(fe_builder.n_blocks());
Expand Down Expand Up @@ -200,7 +200,7 @@ void Stokes<dim>::setup_dofs (const bool &first_run)
DoFTools::make_hanging_node_constraints (*dof_handler,
constraints);

dirichlet_bcs.interpolate_boundary_values(*dof_handler, constraints);
dirichlet_bcs.interpolate_boundary_values(*mapping,*dof_handler, constraints);
constraints.close ();

constraints_dot.clear ();
Expand All @@ -209,7 +209,7 @@ void Stokes<dim>::setup_dofs (const bool &first_run)
DoFTools::make_hanging_node_constraints (*dof_handler,
constraints_dot);

dirichlet_dot.interpolate_boundary_values(*dof_handler, constraints_dot);
dirichlet_dot.interpolate_boundary_values(*mapping,*dof_handler, constraints_dot);
constraints_dot.close ();

jacobian_matrix.clear();
Expand All @@ -219,10 +219,9 @@ void Stokes<dim>::setup_dofs (const bool &first_run)
Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
for (unsigned int c=0; c<dim+1; ++c)
for (unsigned int d=0; d<dim+1; ++d)
if (! ((c==dim) && (d==dim)))
coupling[c][d] = DoFTools::always;
else
coupling[c][d] = DoFTools::none;
coupling[c][d] = DoFTools::always;

coupling[dim][dim] = DoFTools::none;

DoFTools::make_sparsity_pattern (*dof_handler,
coupling,
Expand All @@ -244,9 +243,12 @@ void Stokes<dim>::setup_dofs (const bool &first_run)
for (unsigned int d=0; d<dim+1; ++d)
if (c==d)
prec_coupling[c][d] = DoFTools::always;
else if (c<dim && d<dim)
prec_coupling[c][d] = DoFTools::always;
else
prec_coupling[c][d] = DoFTools::none;


DoFTools::make_sparsity_pattern (*dof_handler,
prec_coupling,
jacobian_preconditioner_matrix_sp,
Expand All @@ -266,8 +268,8 @@ void Stokes<dim>::setup_dofs (const bool &first_run)

if (first_run)
{
VectorTools::interpolate(*dof_handler, initial_solution, solution);
VectorTools::interpolate(*dof_handler, initial_solution_dot, solution_dot);
VectorTools::interpolate(*mapping,*dof_handler, initial_solution, solution);
VectorTools::interpolate(*mapping,*dof_handler, initial_solution_dot, solution_dot);
}

computing_timer.exit_section();
Expand All @@ -284,15 +286,15 @@ void Stokes<dim>::update_constraints (const double &t)
DoFTools::make_hanging_node_constraints (*dof_handler,
constraints);

dirichlet_bcs.interpolate_boundary_values(*dof_handler, constraints);
dirichlet_bcs.interpolate_boundary_values(*mapping,*dof_handler, constraints);

constraints.close ();
constraints_dot.clear();

DoFTools::make_hanging_node_constraints (*dof_handler,
constraints_dot);

dirichlet_dot.interpolate_boundary_values(*dof_handler, constraints_dot);
dirichlet_dot.interpolate_boundary_values(*mapping,*dof_handler, constraints_dot);
constraints_dot.close ();
}

Expand Down Expand Up @@ -327,7 +329,6 @@ void Stokes<dim>::assemble_jacobian_matrix(const double t,

FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
FullMatrix<double> cell_prec (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);

const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
Expand Down Expand Up @@ -501,10 +502,6 @@ int Stokes<dim>::residual (const double t,
const FEValuesExtractors::Vector u (0);
const FEValuesExtractors::Scalar p (dim);

std::vector<Tensor<1,dim> > phi_u(dofs_per_cell);
std::vector<SymmetricTensor<2,dim> > grads_phi_u(dofs_per_cell);
std::vector<double> div_phi_u(dofs_per_cell);

typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler->begin_active(),
endc = dof_handler->end();
Expand Down Expand Up @@ -603,8 +600,7 @@ template <int dim>
void Stokes<dim>::output_step(const double t,
const VEC &solution,
const VEC &solution_dot,
const unsigned int step_number,
const double /* h */ )
const unsigned int step_number)
{
computing_timer.enter_section ("Postprocessing");

Expand Down Expand Up @@ -638,8 +634,6 @@ void Stokes<dim>::output_step(const double t,

template <int dim>
bool Stokes<dim>::solver_should_restart (const double t,
const unsigned int ,
const double ,
VEC &solution,
VEC &solution_dot)
{
Expand Down Expand Up @@ -737,7 +731,6 @@ template <int dim>
int Stokes<dim>::setup_jacobian (const double t,
const VEC &src_yy,
const VEC &src_yp,
const VEC &,
const double alpha)
{
computing_timer.enter_section (" Setup Jacobian");
Expand All @@ -750,12 +743,7 @@ int Stokes<dim>::setup_jacobian (const double t,
}

template <int dim>
int Stokes<dim>::solve_jacobian_system (const double ,
const VEC &,
const VEC &,
const VEC &,
const double ,
const VEC &src,
int Stokes<dim>::solve_jacobian_system (const VEC &src,
VEC &dst) const
{
computing_timer.enter_section (" Solve system");
Expand Down Expand Up @@ -848,7 +836,53 @@ void Stokes<dim>::run ()
constraints.distribute(solution);
constraints_dot.distribute(solution_dot);

dae.start_ode(solution, solution_dot, max_time_iterations);
ida.create_new_vector = [this]() ->shared_ptr<VEC>
{
return this->create_new_vector();
};
ida.residual = [this](const double t,
const VEC &y,
const VEC &y_dot,
VEC &residual) ->int
{
return this->residual(t,y,y_dot,residual);
};

ida.setup_jacobian = [this](const double t,
const VEC &y,
const VEC &y_dot,
const double alpha) ->int
{
return this->setup_jacobian(t,y,y_dot,alpha);
};

ida.solver_should_restart = [this](const double t,
VEC &y,
VEC &y_dot) ->bool
{
return this->solver_should_restart(t,y,y_dot);
};

ida.solve_jacobian_system = [this](const VEC &rhs,
VEC &dst) ->int
{
return this->solve_jacobian_system(rhs,dst);
};

ida.output_step = [this](const double t,
const VEC &y,
const VEC &y_dot,
const unsigned int step_number)
{
this->output_step(t,y,y_dot,step_number);
};

ida.differential_components = [this]() ->VEC &
{
return this->differential_components();
};

ida.solve_dae(solution, solution_dot);

computing_timer.print_summary();
timer_outfile.close();
Expand Down
3 changes: 2 additions & 1 deletion examples/heat_equation/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,8 @@ SET(_main source/main.cc)
# Normally you shouldn't need to change anything below.
############################################################
# Declare all source files the target consists of:
file(GLOB _files source/*cc)
file(GLOB _files source/*cc
include/*.h)

# Don't compile the main file into the library
IF(NOT "${_main}" STREQUAL "")
Expand Down
Loading

0 comments on commit c6d0651

Please sign in to comment.