I am not very familiar with FEM and have limited knowledge about its implementation for the Navier-Stokes equation. Currently, I am trying to understand the principles of FEM and parallel computing, but my learning progress is slow.
For now, my research involves solving 2D/3D N-S equation and temperature equation under a large range of Ra/Re number(10^5<=Ra<=10^10). Consequently, the computational cost increases significantly due to the Kolmogorov length scale. To study a 2D convective heat transfer problem, I tried to adapt the script Rayleigh-Benard-time.edp by frederichecht, it works well but, slower and slower when Ra increases.
My serial script uploaded:R-B_obstacles.edp (7.8 KB)
My supervisor recommends the projection algorithm and parallel computing.
The projection method avoids solving a large system at a time, replaced by 2 or 3 decoupled equations for velocity and a Poisson-like equation for pressure. A implementation NSprojection.edp uses only P1 elements, but some common practices using the pressure space 1 order lower than the velocity space such as P2-P1.
I watched the tutorial Pierre Jolivet introduction to parallel FreeFEM part 1, I noticed that the equation for pressure pb4p(q, w, solver=CG, init=n, eps=epsp) = int2d(Th)(dx(q)*dx(w) + dy(q)*dy(w)) - int2d(Th)((dx(u) + dy(v))*w/dt) could be solved in parallel just as the video did.
To find some power law through a series of simulations under a large range of Ra, I need to ensure my numerical results reflect the properties of the equation derived from true physical world correctly.
So the problems are:
Can NSprojection.edp handle high Reynolds numbers and be easily adapted in parallel?
Hello,
I am not an expert on parallel issues, but I can give you some hints on the other parts.
2. You should not use convect, which is relevant when there is no viscosity. Since you have viscosity you can solve the convection implicitly as
In the projection method it is not mandatory to satisfy an inf-sup condition (thus you can use P1/P1), but you will get a more accurate solution if the inf-sup condition is satisfied (as for P2/P1).
Sorry but I am a bit confused. Why the built-in function convect relevant to inviscid flow?
I have no idea because the convection term u'*grad(u) is without nu.
Also, those serial examples related to N-S equation like NSUzawaCahouetChabart.edp using “convect” and they are indeed viscous flow. Could you please explain this further or if I have some misunderstanding?
Thank you very much for your help,
H. Weng
The fact is that when there is no viscosity the problem is hyperbolic. Thus you should not use the implicit formulation because it could be unstable, and give weak or strong oscillations. You should use convect, which is adapted to the hyperbolic structure.
With viscosity you have elliptic regularization of the gradient of the solution, that stabilizes the implicit resolution of convection (intuitively second-order derivatives dominate the first-order ones involved in the convection). This holds with a moderate Reynolds number Re=L*u/nu. But you are right in this case you can also use convect even if it is of really different nature than the spirit of finite elements (but it is truly relevant if you consider small viscosity i.e. large Reynolds number).
Indeed since there is space discretization the real key value is the discrete Reynolds number Re_h=h*u/nu with h the size of the mesh cell. If Re_h is not too large then the implicit formulation is possible.
Thank you very much for such detailed response. I need some time to digest those different stability properties of PDEs.
So, if the discrete Reynolds number Re_h is not too large, computational stability can be guaranteed, right? From the definition Re_h=h*u/nu, it appears to imply that the size h is proportional to nu/u, when I take Re_h fixed. Before, when simulating some 3D pipe flow I just build a initial 2D mesh, and for different Re=1/nu, and do splitting to the initial mesh to ensure h_local*Re^3/4=const. My problem definition, similar to NSprojection.edp, using convect, but with an extra dimension z, I noticed that sometimes, during the transition to turbulence, the computation blows up when timestep too large.
Is it better, for the stability and accuracy of computing, to use something like adaptmesh to let local mesh size h = const*nu/u? Here, the exponent of Re_h is 1, whereas that of the Kolmogorov scale is 3/4. I would like to know if maintaining the condition of an appropriate Re_h is more strict than that of an appropriate Kolmogorov scale to resolve the details. Further more, what are the criteria for selecting an appropriate mesh size in DNS within FEM?
The condition say Re_h<1 does not mean that you get an accurate solution, it just means that your computation is stable with respect to space discretisation if you use the implicit method for advection.
If you use convect this condition is not required.
If you have other criterion on h to achieve accuracy, you have to respect them additionally.
Concerning time stability, probably you need to respect that nu*dt/h^2 and u*dt/h of order of unity at most.
The mesh size h should be dictacted a priori by your data: it should be the cutoff length under which you do not wish to compute the details of the solution. But as you say it should be small enough in relation with nu in order to resolve your Kolmogorov scale.
I have written a FreeFEM code to compute vorticity formulation of incompressible flow using problem, which is a bit slow and can’t be parallelized.
So I tried to write a matrix form but met some problem.
I found an algorithm to solve this problem(JIAN-GUO LIU AND WEINAN E, 2000"SIMPLE FINITE ELEMENT METHOD IN VORTICITY FORMULATION FOR INCOMPRESSIBLE FLOWS"), which decouples the vorticity \omega and stream function \psi :
The X_{0,h}^k is the subspace of X_h^k with zero boundary value.
To compute an auxiliary term \bar{\omega}^{n+1}, I need to get a submatrix (A_{0,0}A_{0,b}) from the stiffness matrix A to compute right hand side of (I’), and compute the nonlinear term N.
In my code using problem, the explicit nonlinear term (\nabla \varphi,\omega \nabla^{\perp}\psi) is:
grad(v1)'*[-dy(psiold),dx(psiold)]*omegaold
How to write a varf of nonlinear term as (2.6)?
I know that I can use varf Stiffness(u,v)=int2d(Th)(grad(u)'*grad(v)); and set tgv=-10 to build a matrix of Stiffness with zero rows of label 1 , but I have no idea to build A_{0,b}.
I found the method to get list of DoF boundary in Frédéric Hecht’s reply:
If I get the list of DoF on boundary, should I permutate the rows and columns of Stiffness to get A？Or is there a better way to implement such algorithm?
Could you give me some advice?
On the question of extracting a submatrix, you have to make the list of indices corresponding to each block, then you can extract the submatrices. An example is matrix-extract-bdry.edp (1.3 KB)
Dear Prof. Bouchut,
The vorticity formulation gives somewhat abnormal results, which conflicted with previous ones. I don’t have much time left to complete my paper, so I go back to solve the original N-S equation, but in parallel. I am working on a convection problem, unlike Boussinesq approximation, the density is modelled nonlinearly: \rho=\rho_0(1-\beta|T-T_M|^q),where β > 0 is the thermal expansion coefficient and ρ_0 is the density at T = T_M.
Conventionally, q=2, and the dimensionless equation set reads
,other terms are omitted. I found in varf, a term should be either bilinear or linear, so the uT^2*v2 referred to trilinear term, is illegal.
I tried to substitute uT^2*v2 to uT*uTb*v2, uTb is temperature at previous step, and I can run code. If it is a right doing? I also tried to use fully explicit term (uTb-TM)^2, but instability won’t be observed: the system remains stable at a high Rayleigh number, and the growing of time step does not have an effect.
Your two tries uT*uTb*v2 and (uTb-TM)^2 could be possible, numerical stability analysis would be necessary in order to know what is the best. In practice you can just try, see, and choose what you prefer. Next if the system does not behave numerically as you expect (unstable whereas you expect stable, or stable whereas you expect unstable), this is a hard problem of numerical analysis!
If you are NOT working on developing a FEM code in FreeFEM from scratch anf your aim is solve problems using FEM over the mid-range problem, I suggest checking out NEK5000. You can define your problem and run DNS as well as LES for wide-range of Re using high-accuracy and faster computations. It is a well-established software.
NOTE: Just a suggestion, if your aim is to study a problem and not develop code from scratch.
Thanks for the suggestion. I know a bit about Nektar++, but not a skilled user. For the problem I am working on, it is hard for me to handle the complex mesh, and the multi-physical domain in Nektar++. Nek5000 or Nektar should have the competence to solve my problem, but at this stage, I don’t have much time to adapt my FreeFEM++ code to Nektar++.
Also, it takes time to be accustomed to the way Nektar++ solving a problem: to configure a .xml file and run Nektar++ session.
I am currently switching from a IBM based DNS code to NEK5000. It worked well for me. In NEK5000 or Nektar, you don’t have to code the NS equation as it is already there. You just have to tune some parameters. So once you get a hang of it, it is way faster than writing the code from scratch and optimizing it.
