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_pointhttp://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
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
namespace Non_Linear_Poisson
{
using namespace std;
using namespace dealii;
template
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 triangulation;
DoFHandler dof_handler;
FE_Q fe;
ConstraintMatrix hanging_node_constraints;
PETScWrappers::MPI::SparseMatrix system_matrix;
PETScWrappers::MPI::Vectorpresent_solution;
PETScWrappers::MPI::Vectornewton_update;
PETScWrappers::MPI::Vectorsystem_rhs;
};
template
class BoundaryValues : public Function
{
public:
BoundaryValues(): Function () {}
virtual double value(const Point &p, const unsigned int component = 0) const;
};
template
double BoundaryValues::value(const Point &p, const unsigned int /*component*/) const
{
return sin(3*p[0]) * sin(3*p[1]);
}
template
class RightHandSide : public Function
{
public:
RightHandSide(): Function () {}
virtual double value(const Point &p, const unsigned int component = 0) const;
};
template
double RightHandSide::value(const Point &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
NLPoisson::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
NLPoisson::~NLPoisson()
{
dof_handler.clear();
}
template
void NLPoisson::setup_system(const bool initial_step)
{
if (initial_step)
{
pcout<<" inside setup 1"< 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"< 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
void NLPoisson::assemble_system ()
{
pcout<<" inside assembly"< quadrature_formula(3);
FEValues 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 u
