Re: [deal.II] Re: Step-35 Poiseuille flow in a 3D pipe

[Martin Kronbichler]

Dear Jiaqi,

Besides the literature suggested by Daniel, you should also read the paper

J.L. Guermond, P. Minev, J. Shen, An overview of projection methods for 
incompressible flows, Comput Methods Appl Mech Engrg 195, 2006.


Section 3.5 explains this question.

As a more general remark to your question I want to remind you that this 
forum is not primarily for asking questions about methods and their 
implications (go read the relevant papers) but rather more technical 
aspects regarding the implementation of finite element discretizations 
in deal.II. See also here: 
https://groups.google.com/forum/#!topic/dealii/GRZMUTLIm2I


Best,
Martin


On 03.08.2016 04:09, Jiaqi ZHANG wrote:

Hello David,

Thank you for the answer. I know Step-35 doesn't implement velocity 
correction, which


is still confusing me. Without correcting the velocity, it is not 
divergence free, does it make


any sense?

Best,
Jiaqi



2016-08-02 16:37 GMT-04:00 David Wells >:


Hi Jiaqi,

I have to admit that it has been about a year since I last really
looked at step-35 and I do not know the answer off the top of my
head for this. I saw that Daniel responded to this question
elsewhere; was that enough to figure things out?

Thanks,
David Wells


On Friday, July 29, 2016 at 9:19:29 PM UTC-4, Jiaqi ZHANG wrote:

Hello David Wells,

Sorry to bother you, but I have a problem about Step-35. I
have posted the question

a few days ago, and no one answered me, so I have to search in
the mailing list and found

that you have been using it a lot. I was wondering if you
could help me, and the following is

my question:

I read the references listed there and some other papers about
projection methods,

and know that after obtaining the intermediate values for
velocity, say v_n, we need to

correct the velocity by u_n = v_n - grad(2dt/3 phi_n), which
is mentioned in the introduction of Step-35,

but it seems that Step-35 only corrects the pressure, and I
cannot find anything in the code that implements the above
equation.

However, when I correct the velocity at every time step, the
algorithm turns out to be divergent (iterating more than 1000
times at some time step).

I am very confused, did I misunderstand the algorithm or the code?

Thanks in advance.


Best,
Jiaqi



在 2015年3月17日星期二 UTC-4下午8:14:21，David Wells写道：

Hi Joaquin,

I used step-35 a lot, so I may be able to assist you here.

For the velocity inflow, are you sure that what you have
will be zero flow at the edges? The XY slices don't look
quite circular to me, so this may be a problem if you
refine the mesh. Since you are using GMSH, it may be
easiest to use GMSH's extrusion abilities to extrude a
circle along the z axis to solve this.

I believe you handled the pressure correctly (as long as
the domain starts at z = 0 and ends at z = 10).

At least as of a few months ago, step-35 could not run in
3D due to a lack of a compressed sparsity pattern. I fixed
this in my own copy of step-35 (see
http://www.github.com/drwells/step-35a) (yeah, I know I
should submit a patch...). I tweaked the preconditioners
to improve performance, so my version might be useful to
look at.

Another thing you might want to look at is how Abner
implemented labels for the different boundaries. According
to my experiments (see

https://github.com/drwells/dealii-step-35a/blob/master/cylinderextruded.geo):
1 is the no-slip edges
2 is the inflow
3 is the outflow
4 is the cylinder (may not apply for you)

If your mesh file does not label things correctly then
things will get weird. I hope that this is helpful.

Thanks,
David Wells

On Tuesday, March 17, 2015 at 6:08:02 PM UTC-4, Joaquin
wrote:

Hi,

I am trying to run step-35 for the unsteady
incompressible viscous flow in a pipe, for 3D case.
The mesh file I used is attached. The graph of results
is attached. The expected results for any time should
be paraboloid (or parabolic in 2D), but, I can not get
that. The changes I did are:

*1. namespace EquationData:

*template 
double Velocity::value (const Point ,
 const unsigned int) const
{
  if (this->comp == 0)
   

[deal.II] Re: Step-35 Poiseuille flow in a 3D pipe

[David Wells]

Hi Jiaqi,

I have to admit that it has been about a year since I last really looked at 
step-35 and I do not know the answer off the top of my head for this. I saw 
that Daniel responded to this question elsewhere; was that enough to figure 
things out?

Thanks,
David Wells

On Friday, July 29, 2016 at 9:19:29 PM UTC-4, Jiaqi ZHANG wrote:
>
> Hello David Wells,
>
> Sorry to bother you, but I have a problem about Step-35. I have posted the 
> question
>
> a few days ago, and no one answered me, so I have to search in the mailing 
> list and found
>
> that you have been using it a lot. I was wondering if you could help me, 
> and the following is 
>
> my question:
>
> I read the references listed there and some other papers about projection 
> methods, 
>
> and know that after obtaining the intermediate values for velocity, say 
> v_n, we need to 
>
> correct the velocity by u_n = v_n - grad(2dt/3 phi_n), which is mentioned 
> in the introduction of Step-35, 
>
> but it seems that Step-35 only corrects the pressure, and I cannot find 
> anything in the code that implements the above equation.
>
> However, when I correct the velocity at every time step, the algorithm 
> turns out to be divergent (iterating more than 1000 times at some time 
> step). 
>
> I am very confused, did I misunderstand the algorithm or the code?
>
> Thanks in advance.
>
>
> Best,
> Jiaqi 
>
>
>
> 在 2015年3月17日星期二 UTC-4下午8:14:21，David Wells写道：
>>
>> Hi Joaquin,
>>
>> I used step-35 a lot, so I may be able to assist you here.
>>
>> For the velocity inflow, are you sure that what you have will be zero 
>> flow at the edges? The XY slices don't look quite circular to me, so this 
>> may be a problem if you refine the mesh. Since you are using GMSH, it may 
>> be easiest to use GMSH's extrusion abilities to extrude a circle along the 
>> z axis to solve this.
>>
>> I believe you handled the pressure correctly (as long as the domain 
>> starts at z = 0 and ends at z = 10).
>>
>> At least as of a few months ago, step-35 could not run in 3D due to a 
>> lack of a compressed sparsity pattern. I fixed this in my own copy of 
>> step-35 (see http://www.github.com/drwells/step-35a) (yeah, I know I 
>> should submit a patch...). I tweaked the preconditioners to improve 
>> performance, so my version might be useful to look at.
>>
>> Another thing you might want to look at is how Abner implemented labels 
>> for the different boundaries. According to my experiments (see 
>> https://github.com/drwells/dealii-step-35a/blob/master/cylinderextruded.geo
>> ):
>> 1 is the no-slip edges
>> 2 is the inflow
>> 3 is the outflow
>> 4 is the cylinder (may not apply for you)
>>
>> If your mesh file does not label things correctly then things will get 
>> weird. I hope that this is helpful.
>>
>> Thanks,
>> David Wells
>>
>> On Tuesday, March 17, 2015 at 6:08:02 PM UTC-4, Joaquin wrote:
>>>
>>> Hi,
>>>
>>> I am trying to run step-35 for the unsteady incompressible viscous flow 
>>> in a pipe, for 3D case. The mesh file I used is attached. The graph of 
>>> results 
>>> is attached.  The expected results for any time should be paraboloid 
>>> (or parabolic in 2D), but, I can not get that. The changes I did are: 
>>>
>>>
>>>
>>> *1. namespace EquationData: *template 
>>> double Velocity::value (const Point ,
>>>  const unsigned int) const
>>> {
>>>   if (this->comp == 0)
>>> {
>>>   const double Um = 1.5;
>>>   const double R  = 2.;
>>>   return Um*(1 - (p(0)*p(0)+p(1)*p(1))/(2*R*R));
>>> }
>>>   else
>>>   return 0.;
>>> }
>>>
>>> *Through the pipe pressure drop **is linear (**I do not how to make 
>>> changes here):* 
>>>
>>> template 
>>> double Pressure::value (const Point ,
>>>  const unsigned int) const
>>> {
>>>   return 10.-p(2);
>>> }
>>>
>>>  
>>> *2. *
>>>   template 
>>>   void
>>>   NavierStokesProjection::
>>>   create_triangulation_and_dofs (const unsigned int n_refines)
>>>   {
>>> GridIn grid_in;
>>> grid_in.attach_triangulation (triangulation);
>>>
>>> {
>>>   std::string filename = "*cilindro.msh*";
>>>   std::ifstream file (filename.c_str());
>>>   Assert (file, ExcFileNotOpen (filename.c_str()));
>>>   grid_in.read_msh (file);
>>> }
>>>
>>> * 3. *
>>> int main()
>>> {
>>>   try
>>> {
>>>   using namespace dealii;
>>>   using namespace Stepw3;
>>>
>>>   RunTimeParameters::Data_Storage data;
>>>   data.read_data ("parameter-file.prm");
>>>
>>>   deallog.depth_console (data.verbose ? 2 : 0);
>>>
>>>   NavierStokesProjection<3> test (data);
>>>   test.run (data.verbose, data.output_interval);
>>> }
>>> .
>>> .
>>> .
>>>   return 0;
>>> }
>>> 4. 
>>> template 
>>>   void NavierStokesProjection::assemble_vorticity (const bool 
>>> reinit_prec)
>>>   {
>>> Assert (dim == 3, ExcNotImplemented());
>>> if (reinit_prec)
>>>

[deal.II] Re: Step-35 Poiseuille flow in a 3D pipe

[Jiaqi ZHANG]

Hello David Wells,

Sorry to bother you, but I have a problem about Step-35. I have posted the 
question

a few days ago, and no one answered me, so I have to search in the mailing 
list and found

that you have been using it a lot. I was wondering if you could help me, 
and the following is 

my question:

I read the references listed there and some other papers about projection 
methods, 

and know that after obtaining the intermediate values for velocity, say 
v_n, we need to 

correct the velocity by u_n = v_n - grad(2dt/3 phi_n), which is mentioned 
in the introduction of Step-35, 

but it seems that Step-35 only corrects the pressure, and I cannot find 
anything in the code that implements the above equation.

However, when I correct the velocity at every time step, the algorithm 
turns out to be divergent (iterating more than 1000 times at some time 
step). 

I am very confused, did I misunderstand the algorithm or the code?

Thanks in advance.


Best,
Jiaqi 



在 2015年3月17日星期二 UTC-4下午8:14:21，David Wells写道：
>
> Hi Joaquin,
>
> I used step-35 a lot, so I may be able to assist you here.
>
> For the velocity inflow, are you sure that what you have will be zero flow 
> at the edges? The XY slices don't look quite circular to me, so this may be 
> a problem if you refine the mesh. Since you are using GMSH, it may be 
> easiest to use GMSH's extrusion abilities to extrude a circle along the z 
> axis to solve this.
>
> I believe you handled the pressure correctly (as long as the domain starts 
> at z = 0 and ends at z = 10).
>
> At least as of a few months ago, step-35 could not run in 3D due to a lack 
> of a compressed sparsity pattern. I fixed this in my own copy of step-35 
> (see http://www.github.com/drwells/step-35a) (yeah, I know I should 
> submit a patch...). I tweaked the preconditioners to improve performance, 
> so my version might be useful to look at.
>
> Another thing you might want to look at is how Abner implemented labels 
> for the different boundaries. According to my experiments (see 
> https://github.com/drwells/dealii-step-35a/blob/master/cylinderextruded.geo
> ):
> 1 is the no-slip edges
> 2 is the inflow
> 3 is the outflow
> 4 is the cylinder (may not apply for you)
>
> If your mesh file does not label things correctly then things will get 
> weird. I hope that this is helpful.
>
> Thanks,
> David Wells
>
> On Tuesday, March 17, 2015 at 6:08:02 PM UTC-4, Joaquin wrote:
>>
>> Hi,
>>
>> I am trying to run step-35 for the unsteady incompressible viscous flow 
>> in a pipe, for 3D case. The mesh file I used is attached. The graph of 
>> results 
>> is attached.  The expected results for any time should be paraboloid (or 
>> parabolic in 2D), but, I can not get that. The changes I did are: 
>>
>>
>>
>> *1. namespace EquationData: *template 
>> double Velocity::value (const Point ,
>>  const unsigned int) const
>> {
>>   if (this->comp == 0)
>> {
>>   const double Um = 1.5;
>>   const double R  = 2.;
>>   return Um*(1 - (p(0)*p(0)+p(1)*p(1))/(2*R*R));
>> }
>>   else
>>   return 0.;
>> }
>>
>> *Through the pipe pressure drop **is linear (**I do not how to make 
>> changes here):* 
>>
>> template 
>> double Pressure::value (const Point ,
>>  const unsigned int) const
>> {
>>   return 10.-p(2);
>> }
>>
>>  
>> *2. *
>>   template 
>>   void
>>   NavierStokesProjection::
>>   create_triangulation_and_dofs (const unsigned int n_refines)
>>   {
>> GridIn grid_in;
>> grid_in.attach_triangulation (triangulation);
>>
>> {
>>   std::string filename = "*cilindro.msh*";
>>   std::ifstream file (filename.c_str());
>>   Assert (file, ExcFileNotOpen (filename.c_str()));
>>   grid_in.read_msh (file);
>> }
>>
>> * 3. *
>> int main()
>> {
>>   try
>> {
>>   using namespace dealii;
>>   using namespace Stepw3;
>>
>>   RunTimeParameters::Data_Storage data;
>>   data.read_data ("parameter-file.prm");
>>
>>   deallog.depth_console (data.verbose ? 2 : 0);
>>
>>   NavierStokesProjection<3> test (data);
>>   test.run (data.verbose, data.output_interval);
>> }
>> .
>> .
>> .
>>   return 0;
>> }
>> 4. 
>> template 
>>   void NavierStokesProjection::assemble_vorticity (const bool 
>> reinit_prec)
>>   {
>> Assert (dim == 3, ExcNotImplemented());
>> if (reinit_prec)
>>  .
>>  .
>>
>> Into the parameter file: 
>> - Set n_of_refines = 1
>> - Set max_iterations = 5000  # maximal number of iterations that GMRES 
>> must make
>>
>> Please, can someone give me some guidelines to achieve my goal?, I have 
>> one week to present this easy step (for me it is difficult). I have studied 
>> the projection step method, but, the details I can't understand very well 
>> yet. (The boundary conditions are: at r = 0 then v = vmax; r = R then v = 
>> 0, Initial condition: t = 0 then v = 0).