Numerical modeling of Geophysical Flows by Finite Element techniques - PowerPoint PPT Presentation

Numerical modeling of Geophysical Flows by Finite Element techniques with FreeFem++ F. Hecht Laboratoire Jacques-Louis Lions Universit´ e Pierre et Marie Curie Paris, France with O. Pironneau http://www.freefem.org mailto:hecht@ann.jussieu.fr With the support of ANR (French gov.) ANR-07-CIS7-002-01 http://www.freefem.org/ff2a3/ http://www-anr-ci.cea.fr/ Cours ,University of Seville, Spain, March 2010 1

PLAN – Some Fluids equation – Mesh generation – The equation potentiel flow. – Poisson equation with matrix – The Variational formulation. – Variational/Weak form (Matrix and – The FEM in few minutes. vector ) – How to solve this equation in – Time dependent problem FreeFem++ – Stokes variational Problem – Introduction Freefem++ – Navier-Stokes (characteristic me- – some syntaxe thod) – The potentiel flow equation in – Navier-Stokes (Newton method) 2d,3d http://www.freefem.org/ Cours ,University of Seville, Spain, March 2010 2

Basic fluids equations First we suppose the fluids are in a domain Ω ⊂ R d where d = 2 or 3 is the dimension of the physical space. NOTATION : ∂ i u = ∂u x = ( x i ) d i =1 ∈ R d , ∇ u = ( ∂ i u ) d � , a . b = a i b i , i =1 . ∂x i i In lot of the case the fluid are incompressible so if u = ( u i ) d i =1 is the velocity field of the fluid then ∇ . u = 0 = � d i =1 ∂ i u i . If we suppose the fluid are irrotational, we get the existence of a potential field φ such that u = ∇ φ = ( ∂ i φ ) d i =1 . By mixed the both equations we get the potentiel flow equation (also call Poisson equation) : − ∆ φ = −∇ . ∇ φ = 0 (1) where ∆ φ = � d i =1 ∂ 2 ii φ . Cours ,University of Seville, Spain, March 2010 3

Basic fluids equations, BC The boundary condition on the border Γ of Ω of Poisson equation can be φ given (call Dirichlet boundary condition) αφ + ∂ n φ given (call Robin or Fourier boundary condition if α > 0 and Neumann if α = 0 ), where n is the outer unit normal of Γ and ∂ n = n . ∇ . Around a wing Execute Domain-Naca012.edp In a lac Execute Domain-Leman.edp Cours ,University of Seville, Spain, March 2010 4

Poisson equation, weak form Let a domain Ω with a partition of ∂ Ω in Γ d , Γ r . Find φ a solution in such that : αφ + ∂φ − ∆ φ = f in Ω , φ = g d on Γ d , n = g n on Γ r (2) ∂� Remaind Integration by part also call Stokes formula : � � � − Ω ( ∇ . u ) v dΩ = Ω u . ∇ v dΩ − Γ u . n v dΓ Denote V g = { ψ ∈ H 1 (Ω) /ψ | Γ d = g } . The Basic variationnal formulation after multiply by ψ integrated and integrated by part is : find φ ∈ V g d (Ω) , such that ∂φ � � � Ω ∇ φ. ∇ ψ = ∀ ψ ∈ V 0 (Ω) Ω fψ + ∂nψ, (3) Γ Cours ,University of Seville, Spain, March 2010 5

Poisson equation, weak form So the border term in green become ∂φ � � � g n ψ − + ∂nψ = + αφψ Γ Γ r Γ r and the equivalent problem is Find φ ∈ V g d (Ω) , such that � � � � Ω ∇ φ. ∇ ψ + αφψ = Ω fψ + g n ψ, ∀ ψ ∈ V 0 (Ω) (4) Γ r Γ r Cours ,University of Seville, Spain, March 2010 6

Poisson equation, weak full Newann BC Remark, if Γ r = ∅ and α = 0 then φ is defined only through derivative, so φ is defined to a constant,to fixe the constant we add the constraint � Ω φ = 0. Now the system is not square, and if we put ψ = 1 in (4) we get � Ω f + � Γ g n = 0 We can make solve a small perturbated problem with a small positive ε : Find φ ε ∈ H 1 (Ω) � � � ∀ ψ ∈ H 1 (Ω) Ω ∇ φ ε . ∇ ψ + εφ ε ψ dΩ = Ω fψ dΩ+ Γ g n ψ dΓ , (5) If � Ω f + � Γ g n = 0, we can prove : � ε → 0 φ ε = φ, lim and Ω φ ε dΩ = 0 Cours ,University of Seville, Spain, March 2010 7

Finite Element method Construct a finite dimensional space V h to approach the functional space H 1 (Ω). Where h is a parameter which theoretically go to zero (the mesh size). To build V h first, we build a triangular mesh T h of Ω if d = 2 (example Around a wing Execute Domain-Naca012.edp ) or if d = 3 we build a tetrahedral mesh (exemple a Execute Cube.edp ) So a mesh T h is a ”partition” of Ω h in set of element K (triangle in 2d, or resp. tetraedra in 3d), and and set of border element ( segment in 2d or resp. triangle in 3d). We have and K ′ ∩ K = ∅ if K ′ � = K ∈ T h ∪ Ω h = K, K ∈ T h Cours ,University of Seville, Spain, March 2010 8

Finite Element method/ II Let call P k ( K ) is the space of polynomials function of degree ≤ k on K. So now the simple finite element space are : for each k > 0 positive integer we have, V k h = { v ∈ C 0 (Ω h ) / ∀ K ∈ T h , v | K ∈ P k ( K ) } (6) And for the k = 0 the function are constant on K , the function can not be continuous, so the definition is V 0 h = { v ∈ L 2 (Ω h ) / ∀ K ∈ T h , v | K ∈ P 0 ( K ) } (7) Cours ,University of Seville, Spain, March 2010 9

Finite Element method/ III It easy to see that it is finite dimensional vectorial space. We denote w i , i ∈ N k h the basis function of this space V k h , where N k h is the set of degrees of freedom associated to space V k h . We have the following unique decomposition : u h ∈ V k u i w i � h , u h = i ∈N k h and the interpolation operator is defined by ∀ u ∈ C 0 (Ω h ) , I k σ i ( u ) w i � j u = i ∈N k h where σ i are a linear form associated to the degree of freedom which given the value of the degree of freedom. For V 1 h , you can easily prove that : N 1 h is the set of vertices of T h and σ i ( u ) = u ( q i ) (remaind u is a function), where q i is the coordinate of vertex i . Cours ,University of Seville, Spain, March 2010 10

Finite Element method/ IIII To get a approximation of (4) Find φ ∈ V g d (Ω) , such that � � � � Ω ∇ φ. ∇ ψ + ∀ ψ ∈ V 0 (Ω) αφψ = Ω fψ + g n ψ, (8) Γ r Γ r We just replace the space V g by the corresponding finite element space V h g where V h g First defined N 0 h the subset N h with are not the border Γ d , (c-a-d : If the support of σ i is include in Γ d then i is not in N 0 h ) and then we can compute σ i ( g d ) so the definition V g h is : V h g = { ψ h ∈ V h : ∀ i ∈ N h \ N 0 h σ i ( ψ h ) = σ i ( g ) } The discrete problem is : Find φ h ∈ V h g d (Ω) , such that � � � � ∇ φ h . ∇ ψ h + ∀ ψ h ∈ V h 0 (Ω) αφ h ψ h = fψ h + g n ψ h , (9) Ω h Γ he Ω h Γ hr Cours ,University of Seville, Spain, March 2010 11

Introduction to install, freefem++ FreeFem++ is a software to solve numerically partial differential equations R 2 ) and in I R 3 ) with finite elements methods. We used a user (PDE) in I language to set and control the problem. The FreeFem++ language allows for a quick specification of linear PDE’s, with the variational formulation of a linear steady state problem and the user can write they own script to solve no linear problem and time depend problem. You can solve coupled problem or problem with moving domain or eigenvalue problem, do mesh adaptation , compute error indicator, etc ... FreeFem++ is a freeware and this run on Mac, Unix and Window architecture, in parallel with MPI. Fist install freefem++ from http://www.freefem.org/ff++/ftp//FreeFem++-3.8.exe and notepad++ from http://sourceforge.net/project/showfiles.php?group_id=95717&package_id=102072 ‘ Cours ,University of Seville, Spain, March 2010 12

Configuration de Notepad++ Open Notepad++ and Enter F5 In the new window enter the command FreeFem++ "$(FULL_CURRENT_PATH)" Click on Save, and enter FreeFem++ in the box "Name", now choose the short cut keyboard to launch directly FreeFem++ for example alt+shitt+R To add the Color Syntax Compatible to FreeFem++ In Notepad++, In Menu "Parameters"->"Configuration of the Color Syntax". In the list "Language" select C++ Add "edp" in the field "add ext" Select "INSTRUCTION WORD" in the list "Description" and in the field "supplementary key word", cut and past the following list P0 P1 P2 P3 P4 P5 P1dc P2dc P3dc P4dc P5dc RT0 RT1 RT2 RT3 RT4 RT5 macro plot int1d int2d solve movemesh adaptmesh trunc checkmovemesh on func buildmesh square Eigenvalue min max imag exec LinearCG NLCG Newton BFGS LinearGMRES catch try intalledges jump average mean load savemesh convect abs sin cos tan atan asin acos cotan sinh cosh tanh cotanh atanh asinh acosh pow exp log log10 sqrt dx dy endl cout Select "TYPE WORD" in the list "Description" and ... " "supplementary key word", cut and past the following list mesh real fespace varf matrix problem string border complex ifstream ofstream Click on Save& Close. Now nodepad++ is configure. Cours ,University of Seville, Spain, March 2010 13