Re: [deal.II] Parallel Nonlinear Poisson code in DEALII

2019-02-21 Thread m boron 
Dear Jean-Paul,
Thank you very much for the quick response. Yes, it gives an error during
execution. I will also look at the Step-40 tutorial also.

Yours sincerely
Boron

On Thu, Feb 21, 2019 at 6:03 PM Jean-Paul Pelteret 
wrote:

> Dear Boron,
>
> As a new user, I’d like to welcome to the forum and deal.II!
>
> Am I correct that you only see this issue when executing your program in
> parallel? I think that the problem results from the way that you initialise
> the solution vector. That is, this line:
>
> present_solution.reinit(locally_owned_dofs, mpi_communicator);
>
>
> This is perfectly fine in many situations. However, in order to compute
> the solution values and gradients at quadrature points in a cell all
> degrees of freedom belonging to that cell must be accessible from the
> solution vector. Since this vector “ present_solution.reinit" does not
> contain the ghost DoFs, there will eventually be a cell that has its DoFs
> split between processes and not all values will be immediately accessible
> on a process. The way to fix this issue would be to initialise an
> intermediate vector that also holds ghost values, copy the distributed
> solution to the ghosted one, and then proceed with all cell-based
> calculations using the ghosted vector. That would be done something like
> this:
>
>   // At the beginning of assemble_system()
>   PETScWrappers::MPI::Vector present_solution_ghosted;
>   present_solution_ghosted.reinit(locally_owned_dofs,
> locally_relevant_dofs , mpi_communicator);
>   …
>   // On a cell
>fe_values.get_function_values(present_solution_ghosted, old_solution);
>
> A tutorial that shows how to get the locally_relevant_dofs is step-40
> .
>
> I hope that this helps you solve the problem.
> Best,
> Jean-Paul
>
> On 21 Feb 2019, at 11:02, mboron1...@gmail.com wrote:
>
> 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 {
> // 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.
> http://nonlinearpoissonparallel.cc>>
>
>
> --
> 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.
>

-- 
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.

Re: [deal.II] Parallel Nonlinear Poisson code in DEALII

2019-02-21 Thread Jean-Paul Pelteret 
Dear Boron,

As a new user, I’d like to welcome to the forum and deal.II!

Am I correct that you only see this issue when executing your program in 
parallel? I think that the problem results from the way that you initialise the 
solution vector. That is, this line:

> present_solution.reinit(locally_owned_dofs, mpi_communicator);


This is perfectly fine in many situations. However, in order to compute the 
solution values and gradients at quadrature points in a cell all degrees of 
freedom belonging to that cell must be accessible from the solution vector. 
Since this vector “ present_solution.reinit" does not contain the ghost DoFs, 
there will eventually be a cell that has its DoFs split between processes and 
not all values will be immediately accessible on a process. The way to fix this 
issue would be to initialise an intermediate vector that also holds ghost 
values, copy the distributed solution to the ghosted one, and then proceed with 
all cell-based calculations using the ghosted vector. That would be done 
something like this:

  // At the beginning of assemble_system()
  PETScWrappers::MPI::Vector present_solution_ghosted;
  present_solution_ghosted.reinit(locally_owned_dofs, locally_relevant_dofs , 
mpi_communicator);
  …
  // On a cell
   fe_values.get_function_values(present_solution_ghosted, old_solution);  

A tutorial that shows how to get the locally_relevant_dofs is step-40 
.

I hope that this helps you solve the problem.
Best,
Jean-Paul  

> On 21 Feb 2019, at 11:02, mboron1...@gmail.com wrote:
> 
> 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 {
> // 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 
> .
> 

-- 
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.

[deal.II] Parallel Nonlinear Poisson code in DEALII

2019-02-21 Thread mboron1982 
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