Dear all,
I am Boron. I am a new to DEALII. I am currently trying to write a parallel
code in DEALII for solving nonlinear Poisson's equation. The file is also
attahed below. My doubt is "How do we pass history variable while
constructing the cell_matrix?"
A code snippet is (Line No 211-225) :
for (; cell!=endc; ++cell)
{
if (cell->subdomain_id() == this_mpi_process)
{
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);
fe_values.get_function_values(present_solution,
old_solution);
fe_values.get_function_gradients(present_solution,
old_solution_gradients);
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
{
// BUILD ELEMENTAL CELL MATRIX @ EACH GAUSS POINT
}
In the above code snippet, the line
'fe_values.get_function_values(present_solution, old_solution); ' throws
an error. Is there a way to pass only the vectors relevant to the
corresponding subdomains in DEALII?
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
/* PROGRAM TO SOLVE NON-LINEAR POISSON IN PARALLEL USING PETSC*/
// NABLA DOT ((1+U) (NABLA U)) = F
// U =sin(3x) * sin(3y) in (0,1)^2
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/multithread_info.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/utilities.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/petsc_parallel_vector.h>
#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
#include <deal.II/lac/petsc_precondition.h>
#include <deal.II/lac/petsc_solver.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
#include <iostream>
#include <cmath>
#include <ctime>
#include <fstream>
#include <deal.II/numerics/solution_transfer.h>
namespace Non_Linear_Poisson
{
using namespace std;
using namespace dealii;
template <int dim>
class NLPoisson
{
public:
NLPoisson();
~NLPoisson();
void run();
private:
void setup_system(const bool initial_step);
void assemble_system();
unsigned int solve();
void refine_mesh();
void set_boundary_values();
double compute_residual(const double alpha) const;
double determine_step_length() const;
MPI_Comm mpi_communicator;
const unsigned int n_mpi_processes;
const unsigned int this_mpi_process;
ConditionalOStream pcout;
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
FE_Q<dim> fe;
ConstraintMatrix hanging_node_constraints;
PETScWrappers::MPI::SparseMatrix system_matrix;
PETScWrappers::MPI::Vector present_solution;
PETScWrappers::MPI::Vector newton_update;
PETScWrappers::MPI::Vector system_rhs;
};
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
BoundaryValues(): Function<dim> () {}
virtual double value(const Point<dim> &p, const unsigned int component = 0) const;
};
template<int dim>
double BoundaryValues<dim>::value(const Point<dim> &p, const unsigned int /*component*/) const
{
return sin(3*p[0]) * sin(3*p[1]);
}
template <int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide(): Function<dim> () {}
virtual double value(const Point<dim> &p, const unsigned int component = 0) const;
};
template<int dim>
double RightHandSide<dim>::value(const Point<dim> &p, const unsigned int /*component*/) const
{
return 18 * sin(3*p[0]) * sin(3*p[1])* (sin(3*p[0]) * sin(3*p[1]) + 1) -
9*cos(3*p[1])*cos(3*p[1]) * sin(3*p[0])*sin(3*p[0]) -
9*cos(3*p[0])*cos(3*p[0]) * sin(3*p[1])*sin(3*p[1]);
}
template<int dim>
NLPoisson<dim>::NLPoisson()
:
mpi_communicator (MPI_COMM_WORLD),
n_mpi_processes (Utilities::MPI::n_mpi_processes(mpi_communicator)),
this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)),
pcout (std::cout, (this_mpi_process == 0)),
dof_handler(triangulation), fe(1)
{}
template<int dim>
NLPoisson<dim>::~NLPoisson()
{
dof_handler.clear();
}
template <int dim>
void NLPoisson<dim>::setup_system(const bool initial_step)
{
if (initial_step)
{
pcout<<" inside setup 1"<<endl;
GridTools::partition_triangulation (n_mpi_processes, triangulation);
dof_handler.distribute_dofs(fe);
DoFRenumbering::subdomain_wise (dof_handler);
hanging_node_constraints.clear();
DoFTools::make_hanging_node_constraints (dof_handler, hanging_node_constraints);
hanging_node_constraints.clear();
const std::vector<IndexSet> locally_owned_dofs_per_proc = DoFTools::locally_owned_dofs_per_subdomain(dof_handler);
const IndexSet locally_owned_dofs = locally_owned_dofs_per_proc[this_mpi_process];
present_solution.reinit(locally_owned_dofs, mpi_communicator);
}
pcout<<" inside setup 2"<<endl;
DynamicSparsityPattern dsp(dof_handler.n_dofs(),dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp, hanging_node_constraints, false);
const std::vector<IndexSet> locally_owned_dofs_per_proc = DoFTools::locally_owned_dofs_per_subdomain(dof_handler);
const IndexSet locally_owned_dofs = locally_owned_dofs_per_proc[this_mpi_process];
system_matrix.reinit(locally_owned_dofs, locally_owned_dofs, dsp, mpi_communicator);
newton_update.reinit (locally_owned_dofs, mpi_communicator);
system_rhs.reinit (locally_owned_dofs, mpi_communicator);
}
template <int dim>
void NLPoisson<dim>::assemble_system ()
{
pcout<<" inside assembly"<<endl;
QGauss<dim> quadrature_formula(3);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
const RightHandSide<dim> right_hand_side;
std::vector<double> old_solution(n_q_points);
std::vector<Tensor<1, dim>> old_solution_gradients(n_q_points);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
pcout<<" inside assembly1"<<endl;
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end();
pcout<<" inside assembly2"<<endl;
for (; cell!=endc; ++cell)
{
if (cell->subdomain_id() == this_mpi_process)
{
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);
fe_values.get_function_values(present_solution, old_solution);
fe_values.get_function_gradients(present_solution, old_solution_gradients);
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
{
// BUILD ELEMENTAL STIFFNESS MATRIX @ EACH GAUSS POINT
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
cell_matrix(i,j)+= (
fe_values.shape_grad(i, q_point) *
(1+old_solution[q_point]) *
fe_values.shape_grad(j, q_point)
+
fe_values.shape_grad(i, q_point) *
old_solution_gradients[q_point] *
fe_values.shape_value(j, q_point)
) * fe_values.JxW(q_point);
}
}
// BUILD ELEMENTAL FORCE VECTOR @ EACH GAUSS POINT
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
cell_rhs(i) += (
fe_values.shape_value(i,q_point) *
right_hand_side.value(fe_values.quadrature_point (q_point))
-
fe_values.shape_grad(i, q_point) *
(1+old_solution[q_point]) *
old_solution_gradients[q_point]
) * fe_values.JxW(q_point);
}
}
cell->get_dof_indices (local_dof_indices);
hanging_node_constraints.distribute_local_to_global(cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
}
}
pcout<<" inside assembly3"<<endl;
system_matrix.compress(VectorOperation::add);
system_rhs.compress(VectorOperation::add);
pcout<<" inside assembly4"<<endl;
std::map<types::global_dof_index,double> boundary_values;
//VectorTools::interpolate_boundary_values ( dof_handler,
// 0,
// Functions::ZeroFunction<dim>(dim),
// boundary_values);
//pcout<<" inside assembly5"<<endl;
//MatrixTools::apply_boundary_values ( boundary_values,
// system_matrix,
// newton_update,
// system_rhs,
// false);
pcout<<" getting out assembly6"<<endl;
}
template <int dim>
unsigned int NLPoisson<dim>::solve()
{
pcout<<"inside solver"<<endl;
SolverControl solver_control(newton_update.size(), system_rhs.l2_norm()*1e-8);
PETScWrappers::SolverCG solver(solver_control, mpi_communicator);
PETScWrappers::PreconditionBlockJacobi preconditioner(system_matrix);
solver.solve (system_matrix, newton_update, system_rhs, preconditioner);
Vector<double> localized_newton_update (newton_update);
hanging_node_constraints.distribute(localized_newton_update);
newton_update = localized_newton_update;
const double alpha = determine_step_length();
present_solution.add(alpha, newton_update);
pcout<<"done"<<endl;
return solver_control.last_step();
}
template <int dim>
void NLPoisson<dim>::set_boundary_values()
{
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler, 0, BoundaryValues<dim>(), boundary_values);
for(std::map<types::global_dof_index, double>::const_iterator p=boundary_values.begin(); p!=boundary_values.end(); ++p)
{
present_solution(p->first) = p->second;
}
}
template <int dim>
double NLPoisson<dim>::determine_step_length() const
{
return 1;
}
template <int dim>
void NLPoisson<dim>::run()
{
pcout<<"inside run"<<endl;
unsigned int refinement =0;
bool first_step=true;
const double tolerance = 1e-3;
signed int iteration = -1;
GridGenerator::hyper_cube (triangulation);
triangulation.refine_global(1);
double previous_res = 0;
while (first_step || ((previous_res>tolerance) && iteration<0))
{
iteration=iteration+1;
pcout<<"----------------Iteration Number : "<<iteration<<endl;
if (first_step == true)
{
pcout<<"************Inital Mesh************"<<endl;
setup_system(true);
set_boundary_values();
first_step=false;
}
else
{
++refinement;
pcout<<"************Refined Mesh************"<<endl;
}
for (unsigned int inner_iteration=0; inner_iteration<1; ++inner_iteration)
{
pcout<<" iteration: "<<inner_iteration<<endl;
pcout<<"(";
for (unsigned int p=0; p<n_mpi_processes; ++p)
{
pcout<<(p==0 ? ' ' : '+') << (DoFTools::count_dofs_with_subdomain_association (dof_handler, p));
}
pcout << " )" << std::endl;
assemble_system();
previous_res=system_rhs.l2_norm();
}
}
}
}
int main(int argc, char **argv)
{
try
{
using namespace dealii;
using namespace Non_Linear_Poisson;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
Non_Linear_Poisson::NLPoisson<2> prob;
prob.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
