April 2016 - firedrake - imperial.ac.uk

making a vertical slice with u,v,w in one field
by Colin Cotter 25 Apr '16

25 Apr '16

Hi Firedrakers, I'm revisiting the solver approach in the vertical slice model, to try to incorporate it better into dcore. A vertical slice model is obtained from a 3D model, but you assume that all the fields are independent of y, whilst retaining a fully 3D velocity vector. There are two ways to achieve this: 1) Build a mesh that is only one element wide and periodic in the y-direction. 2) Split the velocity space into a part which is tangential to the slice plane (u,w), and a part which is normal to the slice plane (v). In compatible FEM this means making a 2D H(div) space for the tangential part, and a scalar DG space for v. So far we have been doing 2) in a way that requires rather specialised code that is different from what we use in 3D, which makes the code very bloated and hard to test. First of all, is 1) possible in Firedrake? If not, how difficult would it be to make it possible? This would be my preferred solution. If 1) is out, then I'd like to improve our approach to 2) so that it is less intrusive in our code. I'd really like to make a sum of spaces so that we get a 3D vector with the in plane components coming from an H(div) space and the out-of-plane components coming from a DG space. Then, I could call this whole thing my velocity space and e.g. solve mixed velocity pressure problems using it. Is *this* possible? If not, how difficult would it be to make it possible? all the best --cjc -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

Re: [firedrake] Linear problems with zero right hand side
by Andrew McRae 25 Apr '16

25 Apr '16

Hi Daniel, Nice to see you on here. The "F == 0" syntax is usually reserved for nonlinear problems, though it can also be used for linear problems (which I think is what you have?) However, when using this syntax, F should only contain Functions and TestFunctions, not TrialFunctions. The Jacobian (which corresponds to the LHS if F was linear) is then generated automatically using automatic differentiation of F with respect to an appropriate variable. E.g., in the Burgers example, we have F = (inner((u - u_)/timestep, v) + inner(dot(u,nabla_grad(u)), v) + nu*inner(grad(u), grad(v)))*dx Note that u was declared as a Function. Then a few lines later, while (t <= end): solve(F == 0, u) The second argument to the solve call, u, is the variable to solve for (and where the result gets put), and the AD of F is done with respect to this variable. However, when using the "a == L" syntax for linear problems, a must contain TrialFunctions and TestFunctions (and, optionally, Functions), while L must contain TestFunctions (and usually Functions), and so your other variant works. Best, Andrew On 25 April 2016 at 08:59, Daniel Ruprecht <D.Ruprecht(a)leeds.ac.uk> wrote: > Dear Firedrake team, > > I just began to use firedrake and tried to implement a simple heat > equation to get started. I encountered some strange (strange for me that > is) behaviour when using zero as a right hand side in a linear problem. > > Following the nonlinear Burgers' equation example on the website, I > tried to write the diffusion equation u_t = nu*Delta(u) with backward > Euler as > > a = (inner((u - u_)/timestep, v) + nu*inner(grad(u), grad(v))) * dx > > (removing the nonlinear terms essentially) but when solving > > solve(a == 0, temp) > > this generates a ValueError: Provided residual is not a linear form. > > If instead I put the u_ term on the right hand side > > a = (inner(u, v) + timestep * nu * inner(grad(u), grad(v))) * dx > L = inner(u_, v) * dx > solve(a == L, temp) > > everything works fine. Is this working as intended? I could not find a > comment about this on the website, I was confused because the nonlinear > example seemed to suggest the first variant should be the way to do it. > > I attached the script I used to provide a minimum example. This is with > firedrake/2016-03, i.e. the version from March this year and Python 2.7.9. > > Cheers, Daniel > > -- > Dr. Daniel Ruprecht > University Academic Fellow > School of Mechanical Engineering, Room 4.50 > University of Leeds > Leeds LS2 9JT, UK > Email: d.ruprecht(a)leeds.ac.uk > Phone: +44 (0)113 343 22 01 > >

3 2

Linear problems with zero right hand side
by Daniel Ruprecht 25 Apr '16

25 Apr '16

Dear Firedrake team, I just began to use firedrake and tried to implement a simple heat equation to get started. I encountered some strange (strange for me that is) behaviour when using zero as a right hand side in a linear problem. Following the nonlinear Burgers' equation example on the website, I tried to write the diffusion equation u_t = nu*Delta(u) with backward Euler as a = (inner((u - u_)/timestep, v) + nu*inner(grad(u), grad(v))) * dx (removing the nonlinear terms essentially) but when solving solve(a == 0, temp) this generates a ValueError: Provided residual is not a linear form. If instead I put the u_ term on the right hand side a = (inner(u, v) + timestep * nu * inner(grad(u), grad(v))) * dx L = inner(u_, v) * dx solve(a == L, temp) everything works fine. Is this working as intended? I could not find a comment about this on the website, I was confused because the nonlinear example seemed to suggest the first variant should be the way to do it. I attached the script I used to provide a minimum example. This is with firedrake/2016-03, i.e. the version from March this year and Python 2.7.9. Cheers, Daniel -- Dr. Daniel Ruprecht University Academic Fellow School of Mechanical Engineering, Room 4.50 University of Leeds Leeds LS2 9JT, UK Email: d.ruprecht(a)leeds.ac.uk Phone: +44 (0)113 343 22 01

1 0

Symmetry axis boundary condition
by Fryderyk Wilczynski 24 Apr '16

24 Apr '16

Hello again, I'm having some problems trying to implement a symmetry axis boundary condition for a 2d flow (plane poiseuille). My weak form and approach is described in the pdf. When I try to impose no normal flow through the symmetry axis by setting bc2 = DirichletBC(W.sub(0).sub(0), 0, 1) I get the following exception *** Traceback (most recent call last): File "ns_plane_poiseuille.py", line 104, in <module> solve(F_u == 0, u, bcs = [bc0,bc1,bc2],solver_parameters={'pc_type': 'lu', 'pc_factor_mat_solver_package': 'mumps'} ) # File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", line 120, in solve _solve_varproblem(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving.py", line 164, in _solve_varproblem solver.solve() File "<decorator-gen-295>", line 2, in solve File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/pyop2/profiling.py", line 203, in wrapper return f(*args, **kwargs) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/variational_solver.py", line 190, in solve solving_utils.check_snes_convergence(self.snes) File "/home/cserv1_a/apps/install/firedrake/2016-03-el7/firedrake/lib/python2.7/site-packages/firedrake/solving_utils.py", line 62, in check_snes_convergence %s""" % (snes.getIterationNumber(), msg)) RuntimeError: Nonlinear solve failed to converge after 0 nonlinear iterations. Reason: Inner linear solve failed to converge after 0 iterations with reason: unknown reason (petsc4py enum incomplete?) *** If I run the code without that boundary condition, I get unwanted flow through that symmetry axis.... Any advice on setting the symmetry bc would be greatly appreciated Full code: from __future__ import division from firedrake import * import numpy as np class Configuration(object): def __init__(self, **kwargs): for name, value in kwargs.iteritems(): self.__setattr__(name, value) def __setattr__(self, name, value): """Cause setting an unknown attribute to be an error""" if not hasattr(self, name): raise AttributeError("'%s' object has no attribute '%s'" % (type(self).__name__, name)) object.__setattr__(self, name, value) class MyParameters(Configuration): # Model parameters dt = 1.0e-02 # time step Cn=0.03 # Cahn We=1 # Weber Pe=20000 # Peclet Re=1 # Reynolds gravity = False # switch for gravity Fr2=0.1 # Froude number squared beta = 2 # stabilising term rho1 = Constant(1.0) rho2 = Constant(1.0) drho_hat = 0.5*(1 - (rho2/rho1)) alpha = (rho2-rho1)/(rho2+rho1) class phiBC(DirichletBC): def __init__(self, V, value, sid): super(phiBC,self).__init__(V,value,sid) self.nodes = np.arange(V.dof_dset.total_size, dtype=np.int32) myparameters = MyParameters() # Create mesh and define function spaces m=64 mesh = UnitSquareMesh(m, m) P1 = FunctionSpace(mesh, "Lagrange", 1) P2 = VectorFunctionSpace(mesh, "Lagrange", 2) n = FacetNormal(mesh) # Create mixed space W = P2*P1 # Define functions u = Function(W) # u0 = Function(W) # lagged solution # Dirichlet boundary condition for velocity noslip = Constant((0, 0)) inlet = Function(P2).interpolate(Expression(("0.0","0.25*(x[0] + 1)*(x[0] - 1)"))) bc0 = DirichletBC(W.sub(0), inlet, 4) bc1 = DirichletBC(W.sub(0), noslip, 2) bc2 = DirichletBC(W.sub(0).sub(0), 0, 1) # Initial conditions v0, p0 = u0.split() vel_ini = Expression(("0.0", "0.0")) # Initial conditions for velocity v0.interpolate(vel_ini) def F(u, u0, rho, rho0, rho_hat, W,n, myparameters): # access parameters Re=myparameters.Re # Define test functions chi, theta = TestFunctions(W) v, p = split(u) v0, p0 = split(u0) f3 = inner(dot(v0,nabla_grad(v)), chi)*dx \ + (-p*div(chi) )*dx \ + (1/Re)*inner(nabla_grad(v), nabla_grad(chi) )*dx \ + (1/Re)*inner(nabla_grad(v)[j,i], nabla_grad(chi)[i,j])*dx - (2/(3*Re))*div(v)*div(chi)*dx \ - (1/Re)*inner(inner(chi[i], nabla_grad(v)[i,j]), n[j])*(ds(1)+ds(4)+ds(3)) \ f4 = theta*div(v)*dx F = f3 + f4 return F v, p = u0.split() v.rename("Velocity") p.rename("Pressure") outfile = File("./Results/nsch_res.pvd") parameters["matnest"] = False exact = Function(P2).interpolate(inlet) rho1 = myparameters.rho1 rho2 = myparameters.rho2 F_u = F(u,u, rho1, rho1, rho1, W, n, myparameters) t = 0.0 # Solve the system of equations solve(F_u == 0, u, bcs = [bc0,bc1,bc2],solver_parameters={'pc_type': 'lu', 'pc_factor_mat_solver_package': 'mumps'} ) # u0.assign(u) v, p= u0.split() outfile.write(v, p, time=t) print "L2norm = ", sqrt(assemble(inner(v-exact,v-exact)*dx))

1 0

Re: [firedrake] how to set up random initial conditions?
by David Ham 24 Apr '16

24 Apr '16

Hi Francis, You can set random initial conditions by setting the value of the data array on the function. E.g. f.dat.data[:] = np.random.randn(*f.dat.shape) Note that np.random.randn produces normally distributed values. Check out the contents of np.random to find other possibilities. Regards, David On Fri, 22 Apr 2016 at 19:45 Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello, > > I was thinking of taking my Quasi-Geostrophic demo that I recently > submitted and change the initial conditions to be a random field. Then we > can see the evolution of freely evolving 2D turbulence. Unfortunately, I > am not sure how to pick random initial conditions in firedrake. > > I found issue 514 that suggested doing the following: > > import numpy as np > > > perm_array = np.random.randn(Nx, Ny, Nz) > > > I guess this works if you have a regular mesh. But if I have an irregular > mesh then is there a way to define this for the mesh? > > Cheers, Francis > > >

3 4

Re: [firedrake] Pressure outlet boundary condition in plane poiseuille flow
by Colin Cotter 23 Apr '16

23 Apr '16

Oh yeah, I temporarily forgot about viscosity! I'm so used to developing inviscid models... --cjc On 22 April 2016 at 20:56, Kramer, Stephan <s.kramer(a)imperial.ac.uk> wrote: > On 22/04/16 18:49, Colin Cotter wrote: > > OK, are you trying to match an analytical solution? > > > > I think it depends on whether you integrated by parts in the grad term, > > or the divergence. If you integrated by parts in the grad term, then you > > can get the 0 pressure boundary condition naturally i.e. without setting > > anything. > > > > all the best > > --cjc > > To add to that, the natural ("do nothing") boundary condition leads to a > zero stress condition. This is I think also where the difference between > including and excluding the grad-div term comes from. The tangential > component of stress is different (dv/dx vs. dv/dx+du/dy) and is not > sufficiently constrained by the boundary conditions. If you include the > grad-div term the tangential stress is non-zero and so you either have > to explicitly add that as a boundary term to weakly enforce the correct > value of dv/dx+du/dy, or you apply a strong dirichlet on the tangential > component of velocity. Without the grad-div term the "do-nothing" bc for > the tangential stress automatically enforces dv/dx=0. > > Cheers > Stephan > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

3 3

Re: [firedrake] Pressure outlet boundary condition in plane poiseuille flow
by Colin Cotter 22 Apr '16

22 Apr '16

OK, are you trying to match an analytical solution? I think it depends on whether you integrated by parts in the grad term, or the divergence. If you integrated by parts in the grad term, then you can get the 0 pressure boundary condition naturally i.e. without setting anything. all the best --cjc On 22 April 2016 at 14:58, Fryderyk Wilczynski <scfw(a)leeds.ac.uk> wrote: > I am using P2 for velocity and P1 for pressure, where > > > P1 = FunctionSpace(mesh, "Lagrange", 1) > P2 = VectorFunctionSpace(mesh, "Lagrange", 2) > > thanks > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > on behalf of Colin Cotter <colin.cotter(a)imperial.ac.uk> > *Sent:* 22 April 2016 14:31 > *To:* firedrake > *Subject:* Re: [firedrake] Pressure outlet boundary condition in plane > poiseuille flow > > What kind of finite element spaces are you using? > > all the best > --cjc > > On 22 April 2016 at 12:51, Fryderyk Wilczynski <scfw(a)leeds.ac.uk> wrote: > >> Hi all, >> >> >> I trying to test 2d incompressible navier-stokes code, using plane >> poiseuille flow test case . >> >> I am solving steady state system on a square domain. No-slip on top and >> bottom walls. Velocity inlet on the left where I specify Dirichlet >> parabolic profile ("x[1]*(1-x[1])", "0.0"). On the right boundary I specify >> zero constant pressure. >> >> noslip = Constant((0, 0)) >> inlet = Function(P2).interpolate(Expression(("x[1]*(1-x[1])", "0.0"))) >> bc0 = DirichletBC(W.sub(0), inlet, 1) >> bc1 = DirichletBC(W.sub(1), 0, 2) # zero pressure outlet >> bc2 = DirichletBC(W.sub(0), noslip, 3) >> bc3 = DirichletBC(W.sub(0), noslip, 4) >> >> >> I expect the exact to be v = ( y (1 - y), 0) - (unidirectional flow). >> >> >> However, I think I must have set that outlet condition up incorrectly as >> I get vertical flow near the outlet. Is this a common issue? or a sign that >> there's some bug in my formulation? >> >> >> Thanks, >> >> Fryderyk >> >> >> >> >> > > > -- > http://www.imperial.ac.uk/people/colin.cotter > > www.cambridge.org/9781107663916 > > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1

how to set up random initial conditions?
by Francis Poulin 22 Apr '16

22 Apr '16

Hello, I was thinking of taking my Quasi-Geostrophic demo that I recently submitted and change the initial conditions to be a random field. Then we can see the evolution of freely evolving 2D turbulence. Unfortunately, I am not sure how to pick random initial conditions in firedrake. I found issue 514 that suggested doing the following: import numpy as np perm_array = np.random.randn(Nx, Ny, Nz) I guess this works if you have a regular mesh. But if I have an irregular mesh then is there a way to define this for the mesh? Cheers, Francis

1 0

Re: [firedrake] Pressure outlet boundary condition in plane poiseuille flow
by Colin Cotter 22 Apr '16

22 Apr '16

What kind of finite element spaces are you using? all the best --cjc On 22 April 2016 at 12:51, Fryderyk Wilczynski <scfw(a)leeds.ac.uk> wrote: > Hi all, > > > I trying to test 2d incompressible navier-stokes code, using plane > poiseuille flow test case . > > I am solving steady state system on a square domain. No-slip on top and > bottom walls. Velocity inlet on the left where I specify Dirichlet > parabolic profile ("x[1]*(1-x[1])", "0.0"). On the right boundary I specify > zero constant pressure. > > noslip = Constant((0, 0)) > inlet = Function(P2).interpolate(Expression(("x[1]*(1-x[1])", "0.0"))) > bc0 = DirichletBC(W.sub(0), inlet, 1) > bc1 = DirichletBC(W.sub(1), 0, 2) # zero pressure outlet > bc2 = DirichletBC(W.sub(0), noslip, 3) > bc3 = DirichletBC(W.sub(0), noslip, 4) > > > I expect the exact to be v = ( y (1 - y), 0) - (unidirectional flow). > > > However, I think I must have set that outlet condition up incorrectly as I > get vertical flow near the outlet. Is this a common issue? or a sign that > there's some bug in my formulation? > > > Thanks, > > Fryderyk > > > > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 2

Pressure outlet boundary condition in plane poiseuille flow
by Fryderyk Wilczynski 22 Apr '16

22 Apr '16

Hi all, I trying to test 2d incompressible navier-stokes code, using plane poiseuille flow test case . I am solving steady state system on a square domain. No-slip on top and bottom walls. Velocity inlet on the left where I specify Dirichlet parabolic profile ("x[1]*(1-x[1])", "0.0"). On the right boundary I specify zero constant pressure. noslip = Constant((0, 0)) inlet = Function(P2).interpolate(Expression(("x[1]*(1-x[1])", "0.0"))) bc0 = DirichletBC(W.sub(0), inlet, 1) bc1 = DirichletBC(W.sub(1), 0, 2) # zero pressure outlet bc2 = DirichletBC(W.sub(0), noslip, 3) bc3 = DirichletBC(W.sub(0), noslip, 4) I expect the exact to be v = ( y (1 - y), 0) - (unidirectional flow). However, I think I must have set that outlet condition up incorrectly as I get vertical flow near the outlet. Is this a common issue? or a sign that there's some bug in my formulation? Thanks, Fryderyk

1 0