# HG changeset patch # User gedeone-octave # Date 1384112024 0 # Node ID 945c19831a16b28027a95839edfd8eafce5e7afa # Parent 0c748179f6d495905d4835f8e206bb63ceb5aa6c Unsteady NS with L2 penalization. diff -r 0c748179f6d4 -r 945c19831a16 devel/example/Ficticious_Domain/Unsteady/L2_penalization/PressureUpdate.ufl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/devel/example/Ficticious_Domain/Unsteady/L2_penalization/PressureUpdate.ufl Sun Nov 10 19:33:44 2013 +0000 @@ -0,0 +1,30 @@ +## Copyright (C) 2013 Marco Vassallo +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see . + + +V = VectorElement("CG", triangle, 2) +Q = FiniteElement("CG", triangle, 1) + +# Define trial and test functions +p = TrialFunction(Q) +q = TestFunction(Q) + +# Define coefficients +k = Constant(triangle) +u1 = Coefficient(V) + +# Define bilinear and linear forms +a = inner(grad(p), grad(q))*dx +L = -(1/k)*div(u1)*q*dx diff -r 0c748179f6d4 -r 945c19831a16 devel/example/Ficticious_Domain/Unsteady/L2_penalization/TentativeVelocity.ufl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/devel/example/Ficticious_Domain/Unsteady/L2_penalization/TentativeVelocity.ufl Sun Nov 10 19:33:44 2013 +0000 @@ -0,0 +1,37 @@ +## Copyright (C) 2013 Marco Vassallo +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see . + + +V = VectorElement("CG", triangle, 2) +Q = FiniteElement("CG", triangle, 1) +P0 = FiniteElement("DG", triangle, 0) +# Define trial and test functions +u = TrialFunction(V) +v = TestFunction(V) + +# Define coefficients +k = Constant(triangle) +u0 = Coefficient(V) +p0 = Coefficient(Q) +f = Coefficient(V) +sig = Coefficient(P0) +nu = Coefficient(P0) + +# Define bilinear and linear forms +eq = (1/k)*inner(u - u0, v)*dx + inner(grad(u)*u0, v)*dx + \ + nu*inner(grad(u), grad(v))*dx + sig * (inner (u, v))*dx \ + - inner(f, v)*dx - div(v)*p0*dx +a = lhs(eq) +L = rhs(eq) diff -r 0c748179f6d4 -r 945c19831a16 devel/example/Ficticious_Domain/Unsteady/L2_penalization/UNS_with_L2_penalization.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/devel/example/Ficticious_Domain/Unsteady/L2_penalization/UNS_with_L2_penalization.m Sun Nov 10 19:33:44 2013 +0000 @@ -0,0 +1,128 @@ +## Copyright (C) 2013 Marco Vassallo +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see . + +# In this example we solve the 2D unsteady NS equation for a flow around a square +# cylinder in a channel. In this example, the no-slip BC are computed using +# a L2 penalization technique. +# We use the differential Chorin-Temam algorithm for each time step. +# The solution shows the typical von-karman street pattern after a transition +# period. + +pkg load fem-fenics msh + +x = linspace (0, 10, 80*4+1); +y = linspace (-2, 2, 80+1); +msho = msh2m_structured_mesh (x, y, 1, 1:4); +mesh = Mesh (msho); +plot (mesh); +import_ufl_Problem ("TentativeVelocity"); +import_ufl_Problem ("VelocityUpdate"); +import_ufl_Problem ("PressureUpdate"); + +TH1 = FunctionSpace ('TentativeVelocity', mesh); +TH2 = FunctionSpace ('PressureUpdate', mesh); + +bc1 = DirichletBC (TH1, @(x, y, z, n) [0, 0], [1, 3]); +bc2 = DirichletBC (TH1, @(x, y) [1.5*4*(2 - y)*(2 + y)/16 , 0], 4); +bc3 = DirichletBC (TH2, @(x, y) 0, 2); +bcu = {bc1, bc2}; +bcp = bc3; +## Set parameter values +dt = 0.1; +T = 50.; +u0 = Expression ('u0', @(x, y) [0; 0]); +p0 = Expression ('p0', @(x, y) 0); +k = Constant ('k', dt); +f = Constant ('f', [0; 0]); + +nu_1 = 5.e-3; +nu_0 = 5.e-3; +r = 0.25; +#ficticious domain + +## cylinder +# dom = @(x,y) (((x - 3).^2 + (y).^2) < (r).^2); + +# square +dom = @(x, y) (x <= 3+r)*(x >= 3-r)*(y >= -r)*(y <= r); +nu = Expression ('nu', @(x, y) dom(x, y) * nu_0 + ... + (1 - dom(x, y)) * nu_1); + +sig_1 = 0; +sig_0 = 1e8; +sig = Expression ('sig',@(x, y) dom(x, y) * sig_0 + ... + (1 - dom(x, y)) * sig_1); + + +a2 = BilinearForm ('PressureUpdate', TH2, TH2); +A2 = assemble (a2, bcp); + +a3 = BilinearForm ('VelocityUpdate', TH1, TH1); +A3 = assemble(a3); + +## Time-stepping +t = dt; i = 0; +pold = zeros (26001, 1); +while t < T + + # Compute tentative velocity step + a1 = BilinearForm ('TentativeVelocity', TH1, TH1, k, u0, sig, nu); + L1 = LinearForm ('TentativeVelocity', TH1, f, p0, k, u0); + [A1, b1] = assemble_system (a1, L1, bcu{:}); + utmp = A1 \ b1; + u1 = Function ('u1', TH1, utmp); + + # Pressure correction + L2 = LinearForm ('PressureUpdate', TH2, k, u1); + b2 = assemble (L2, bcp); + ptmp = A2 \ b2; + p1 = Function ('p1', TH2, ptmp); + + # Velocity correction + L3 = LinearForm ('VelocityUpdate', TH1, k, u1, p1); + b3 = assemble (L3); + ut = A3 \ b3; + u1 = Function ('u0', TH1, ut); + + #Update p + p = ptmp + pold; + p0 = Function ('p0', TH2, p); + + # Save to file + save (p0, sprintf ("p_%3.3d", ++i)); + delete (sprintf ("p_%3.3d.pvd", i)); + save (u1, sprintf ("u_%3.3d", i)); + delete (sprintf ("u_%3.3d.pvd", i)); + + # Move to next time step + pold = p; + u0 = u1; + t += dt; + +end + +plot (u0); + +norm_err = 0; +for i = 1:size(msho.p, 2) + x = msho.p (1, i); + y = msho.p (2, i); + if (dom (x, y) == true) + err_L2 = feval (u0, [x, y]); + norm_err += err_L2(1).^2 + err_L2(2).^2; + endif +endfor + +error_on_bc = sqrt (norm_err) diff -r 0c748179f6d4 -r 945c19831a16 devel/example/Ficticious_Domain/Unsteady/L2_penalization/VelocityUpdate.ufl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/devel/example/Ficticious_Domain/Unsteady/L2_penalization/VelocityUpdate.ufl Sun Nov 10 19:33:44 2013 +0000 @@ -0,0 +1,31 @@ +## Copyright (C) 2013 Marco Vassallo +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see . + +# Define function spaces (P2-P1) +V = VectorElement("CG", triangle, 2) +Q = FiniteElement("CG", triangle, 1) + +# Define trial and test functions +u = TrialFunction(V) +v = TestFunction(V) + +# Define coefficients +k = Constant(triangle) +u1 = Coefficient(V) +p1 = Coefficient(Q) + +# Define bilinear and linear forms +a = inner(u, v)*dx +L = inner(u1, v)*dx - k*inner(grad(p1), v)*dx