April 2016 - firedrake - imperial.ac.uk

Re: [firedrake] parallelization
by Andrew McRae 26 Apr '16

26 Apr '16

Dear Francis, What is job.py? If it is a "normal" firedrake script then approach 1) will not be running in parallel. If it is a python script that executes "mpirun -n 16 python simulation.py", on the other hand... Best, Andrew On 25 April 2016 at 18:59, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello, > > So far I have only been running firedrake runs using > > python job.py > > and I thought it was in serial. > > I just tried a run with a higher resolution on my dell that has 2x8 cores > and it seems to be using nearly 32,000% of my cpu. I am very happy to know > that it is parallelizing automatically, it makes my life a lot easier but I > have two questions. > > 1) What kind of parallelization is it using? From the high percentage it > seems like it's threaded but not sure how to check. > > 2) Is there any advantage to running > > mpirun -np 16 python job.py > > versus doing the more simple execution above? > > If not I will continue doing the simple approach. > > Cheers, Francis > > > > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 > >

3 6

Re: [firedrake] Laplace equation on subdomain
by David Ham 26 Apr '16

26 Apr '16

Hi Tomasz, This is a subtlety of UFL. *All* arguments and coefficients in an interior facet integral have to be "restricted". This even applies to ones which are continuous and for which the "+" and "-" restrictions are therefore the same. You can just arbitrarily decide to restrict the continuous things to the plus side, or take avg() of the lack of symmetry offends your sense of aesthetics. Regards, David On Tue, 26 Apr 2016 at 15:46 Tomasz Salwa [RPG] <mmtjs(a)leeds.ac.uk> wrote: > Hi, > > > Previous problem continued. > > I have a working solver for separate subsystems in mixed function space > now. I'd like to introduce the coupling terms at the interface, but they > cause a UFL exception "Discontinuous type Argument must be restricted." > > Full code is here: https://bitbucket.org/mmtjs/lin_coupled_fd_2d > > > <https://bitbucket.org/mmtjs/lin_coupled_fd_2d> > > The relevant part from lib/coupled_system.py is: > > ... > > # define the normal vector at the interface > > self.n = fd.FacetNormal(self.mesh) > self.n_int = (self.B.I("+")*self.W.I("-") - > self.B.I("-")*self.W.I("+"))*self.n("+") # btw how can write this to file > to verify? > > ... > > # the second line causes exception (and later similar terms too) > > a_U = fd.dot( self.trial_s * self.B.I, self.v_s * self.B.I) * > fd.dx \ > + fd.dot( self.v_s, self.n_int ) * self.trial_f * fd.dS > > > Regards, > > Tomasz > > > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > on behalf of Tomasz Salwa [RPG] <mmtjs(a)leeds.ac.uk> > *Sent:* 21 April 2016 11:04 > > *To:* firedrake > *Subject:* Re: [firedrake] Laplace equation on subdomain > > > Cool, works now, dynamic case too. I'll do the same for the structural > part. > > > Cheers, > > Tomasz > > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > on behalf of David Ham <David.Ham(a)imperial.ac.uk> > *Sent:* 21 April 2016 09:38 > *To:* firedrake > *Subject:* Re: [firedrake] Laplace equation on subdomain > Hi Tomasz, > > See the new pull request. In essence, the projection is not safe in the > indicator function case. > > Regards, > > David > > On Wed, 20 Apr 2016 at 15:53 Tomasz Salwa [RPG] <mmtjs(a)leeds.ac.uk> wrote: > >> Hi David, >> >> >> Thank you, now it converges properly. >> >> However, I noticed that the error depends on the mesh resolution beyond >> the fluid subdomain, which I think should not be the case (at least when I >> refine it in x direction). Perhaps the solution is ok, but the source is in >> the way I impose the analytical one, or compute the error. Can you diagnose >> this? >> >> Bitbucket's repo is updated. >> >> >> Regards, >> >> Tomasz >> >> >> ------------------------------ >> *From:* firedrake-bounces(a)imperial.ac.uk < >> firedrake-bounces(a)imperial.ac.uk> on behalf of David Ham < >> David.Ham(a)imperial.ac.uk> >> *Sent:* 20 April 2016 10:05 >> >> *To:* firedrake >> *Subject:* Re: [firedrake] Laplace equation on subdomain >> Hi Tomasz, >> >> I have just put in a pull request fixing your problem. In essence, you >> were only applying half of the indicator function fix. In order to >> correctly use indicator functions to solve problems over a subdomain, it is >> necessary to *both* use the custom boundary condition to force the solution >> in the rest of the domain to zero *and* multiply the integrands in your >> forms by the indicator function. You were only doing the former and not the >> latter. Making this fix causes the error to drop to 10^-5. I haven't done a >> convergence test to check that the error converges at the right rate, but >> you should really do so. >> >> Regards, >> >> David >> >> On Tue, 19 Apr 2016 at 11:52 Tomasz Salwa [RPG] <mmtjs(a)leeds.ac.uk> >> wrote: >> >>> Hi David, >>> >>> >>> Here's the code: >>> >>> https://bitbucket.org/mmtjs/laplace_extended_domain >>> >>> >>> It's roughly the same as before, just with a gmsh-generated mesh. >>> >>> >>> Regards, >>> >>> Tomasz >>> >>> >>> ------------------------------ >>> *From:* firedrake-bounces(a)imperial.ac.uk < >>> firedrake-bounces(a)imperial.ac.uk> on behalf of David Ham < >>> David.Ham(a)imperial.ac.uk> >>> *Sent:* 19 April 2016 09:54 >>> *To:* firedrake >>> >>> *Subject:* Re: [firedrake] Laplace equation on subdomain >>> Hi Tomasz, >>> >>> Can you please post the code you are using for this (preferably by >>> pointing us to a git repo we can clone). I can't diagnose what is going >>> wrong just by looking at the visualisation, I'm afraid. >>> >>> Regards, >>> >>> David >>> >>> On Mon, 18 Apr 2016 at 12:12 Tomasz Salwa [RPG] <mmtjs(a)leeds.ac.uk> >>> wrote: >>> >>>> Hi, >>>> >>>> I plotted the difference between analytical and numerical solutions >>>> (see attached). It doesn't converge to the right solution with mesh >>>> refinement. Is it possible that the Neumann BC is not applied at the >>>> interface during the computation, although the extended Dirichlet BCs >>>> enforce resulting zero solution beyond the fluid subdomain? >>>> >>>> Cheers, >>>> Tomasz >>>> >>>> ________________________________________ >>>> From: firedrake-bounces(a)imperial.ac.uk < >>>> firedrake-bounces(a)imperial.ac.uk> on behalf of Lawrence Mitchell < >>>> lawrence.mitchell(a)imperial.ac.uk> >>>> Sent: 15 April 2016 16:00 >>>> To: firedrake(a)imperial.ac.uk >>>> Subject: Re: [firedrake] Laplace equation on subdomain >>>> >>>> > On 15 Apr 2016, at 15:54, Tomasz Salwa [RPG] <mmtjs(a)leeds.ac.uk> >>>> wrote: >>>> > >>>> > Hi David, >>>> > >>>> > It works partially, i.e. the nodes beyond the fluid subdomain are >>>> "turned off" now, but it looks like the interface is "smeared". I suspect >>>> that the assembled matrix's entries involve also integration over elements >>>> one layer beyond the subdomain boundary. Is it possible to exclude them? >>>> >>>> If you have neumann conditions on your subdomain boundary you've only >>>> constrained the gradient of the solution. Therefore on the boundary nodes >>>> you may have non-zero solution values. >>>> >>>> Because the domain you visualise in paraview extends outside the >>>> subdomain you'll have (in 1D) something like: >>>> >>>> >>>> ........ >>>> \............ >>>> ^ boundary node >>>> >>>> Note how the cell at the boundary has a non-zero coefficient value are >>>> one end, but a zero value at the other (assuming piecewise liners). >>>> Paraview will linearly interpolate this basis function when visualising, >>>> which may be the artefact you're seeing. >>>> >>>> Can you construct an analytic solution and determine that you converge >>>> at the correct rate as you refine the mesh? This is a safer test than >>>> eyeballing the solution. >>>> >>>> Lawrence >>>> >>>

1 0

Re: [firedrake] Laplace equation on subdomain
by David Ham 26 Apr '16

26 Apr '16

Hi Tomasz, See the new pull request. In essence, the projection is not safe in the indicator function case. Regards, David On Wed, 20 Apr 2016 at 15:53 Tomasz Salwa [RPG] <mmtjs(a)leeds.ac.uk> wrote: > Hi David, > > > Thank you, now it converges properly. > > However, I noticed that the error depends on the mesh resolution beyond > the fluid subdomain, which I think should not be the case (at least when I > refine it in x direction). Perhaps the solution is ok, but the source is in > the way I impose the analytical one, or compute the error. Can you diagnose > this? > > Bitbucket's repo is updated. > > > Regards, > > Tomasz > > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > on behalf of David Ham <David.Ham(a)imperial.ac.uk> > *Sent:* 20 April 2016 10:05 > > *To:* firedrake > *Subject:* Re: [firedrake] Laplace equation on subdomain > Hi Tomasz, > > I have just put in a pull request fixing your problem. In essence, you > were only applying half of the indicator function fix. In order to > correctly use indicator functions to solve problems over a subdomain, it is > necessary to *both* use the custom boundary condition to force the solution > in the rest of the domain to zero *and* multiply the integrands in your > forms by the indicator function. You were only doing the former and not the > latter. Making this fix causes the error to drop to 10^-5. I haven't done a > convergence test to check that the error converges at the right rate, but > you should really do so. > > Regards, > > David > > On Tue, 19 Apr 2016 at 11:52 Tomasz Salwa [RPG] <mmtjs(a)leeds.ac.uk> wrote: > >> Hi David, >> >> >> Here's the code: >> >> https://bitbucket.org/mmtjs/laplace_extended_domain >> >> >> It's roughly the same as before, just with a gmsh-generated mesh. >> >> >> Regards, >> >> Tomasz >> >> >> ------------------------------ >> *From:* firedrake-bounces(a)imperial.ac.uk < >> firedrake-bounces(a)imperial.ac.uk> on behalf of David Ham < >> David.Ham(a)imperial.ac.uk> >> *Sent:* 19 April 2016 09:54 >> *To:* firedrake >> >> *Subject:* Re: [firedrake] Laplace equation on subdomain >> Hi Tomasz, >> >> Can you please post the code you are using for this (preferably by >> pointing us to a git repo we can clone). I can't diagnose what is going >> wrong just by looking at the visualisation, I'm afraid. >> >> Regards, >> >> David >> >> On Mon, 18 Apr 2016 at 12:12 Tomasz Salwa [RPG] <mmtjs(a)leeds.ac.uk> >> wrote: >> >>> Hi, >>> >>> I plotted the difference between analytical and numerical solutions (see >>> attached). It doesn't converge to the right solution with mesh refinement. >>> Is it possible that the Neumann BC is not applied at the interface during >>> the computation, although the extended Dirichlet BCs enforce resulting zero >>> solution beyond the fluid subdomain? >>> >>> Cheers, >>> Tomasz >>> >>> ________________________________________ >>> From: firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> >>> on behalf of Lawrence Mitchell <lawrence.mitchell(a)imperial.ac.uk> >>> Sent: 15 April 2016 16:00 >>> To: firedrake(a)imperial.ac.uk >>> Subject: Re: [firedrake] Laplace equation on subdomain >>> >>> > On 15 Apr 2016, at 15:54, Tomasz Salwa [RPG] <mmtjs(a)leeds.ac.uk> >>> wrote: >>> > >>> > Hi David, >>> > >>> > It works partially, i.e. the nodes beyond the fluid subdomain are >>> "turned off" now, but it looks like the interface is "smeared". I suspect >>> that the assembled matrix's entries involve also integration over elements >>> one layer beyond the subdomain boundary. Is it possible to exclude them? >>> >>> If you have neumann conditions on your subdomain boundary you've only >>> constrained the gradient of the solution. Therefore on the boundary nodes >>> you may have non-zero solution values. >>> >>> Because the domain you visualise in paraview extends outside the >>> subdomain you'll have (in 1D) something like: >>> >>> >>> ........ >>> \............ >>> ^ boundary node >>> >>> Note how the cell at the boundary has a non-zero coefficient value are >>> one end, but a zero value at the other (assuming piecewise liners). >>> Paraview will linearly interpolate this basis function when visualising, >>> which may be the artefact you're seeing. >>> >>> Can you construct an analytic solution and determine that you converge >>> at the correct rate as you refine the mesh? This is a safer test than >>> eyeballing the solution. >>> >>> Lawrence >>> >>

2 2

Re: [firedrake] parallelization
by Fabio Luporini 26 Apr '16

26 Apr '16

Hi Francis, We do not support threads, just processes via MPI. What you see is odd, at least to me. So, you're saying that if you run htop, or top, from another shell, you see 1 instance of "python firedrake" taking up 3200% of the CPU (this includes hyperthreads, apparently), arent' you ? What happens if instead you try to run with mpirun ? Can you paste your program somewhere ? -- Fabio 2016-04-26 11:48 GMT+01:00 Francis Poulin <fpoulin(a)uwaterloo.ca>: > Dear Andrew, > > Thanks for the response. > > Sorry. When I said job.py I meant any "normal" firedrake script. When I > run my job on my mac using just python it does run in serial and I have > tested that I can use mpirun to run it in parallel. This is what I > expected. It is on my ubuntu machine that it seems to take up 32 cores > even though I just used python, not mpirun. I wonder whether it might be > using threads instead of cores? > > I have not tried mpirun on this machine yet but I can easily do that. I > am curious to better learn how efficient the parallelization is so I plan > to do tests for different numbers of cores. > > Best regards, > Francis > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] > on behalf of Andrew McRae [A.T.T.McRae(a)bath.ac.uk] > *Sent:* Tuesday, April 26, 2016 6:25 AM > *To:* firedrake(a)imperial.ac.uk > *Subject:* Re: [firedrake] parallelization > > Dear Francis, > > What is job.py? If it is a "normal" firedrake script then approach 1) > will not be running in parallel. If it is a python script that executes > "mpirun -n 16 python simulation.py", on the other hand... > > Best, > Andrew > > On 25 April 2016 at 18:59, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > >> Hello, >> >> So far I have only been running firedrake runs using >> >> python job.py >> >> and I thought it was in serial. >> >> I just tried a run with a higher resolution on my dell that has 2x8 cores >> and it seems to be using nearly 32,000% of my cpu. I am very happy to know >> that it is parallelizing automatically, it makes my life a lot easier but I >> have two questions. >> >> 1) What kind of parallelization is it using? From the high percentage it >> seems like it's threaded but not sure how to check. >> >> 2) Is there any advantage to running >> >> mpirun -np 16 python job.py >> >> versus doing the more simple execution above? >> >> If not I will continue doing the simple approach. >> >> Cheers, Francis >> >> >> >> >> ------------------ >> Francis Poulin >> Associate Professor >> Department of Applied Mathematics >> University of Waterloo >> >> email: fpoulin(a)uwaterloo.ca >> Web: https://uwaterloo.ca/poulin-research-group/ >> Telephone: +1 519 888 4567 x32637 >> >> >

1 0

Saving\loading the solution
by George Ovchinnikov 26 Apr '16

26 Apr '16

Dear all, I'm saving a scalar field as vtu file. I can view is as paraview "function_3" in "Cell\Point Array Status window", but when I try to read it with the following script: import numpy from vtk import * from vtk.util import numpy_support as VN reader = vtkXMLPUnstructuredGridReader() reader.SetFileName('w_0.vtu') # file is attached reader.Update() from pprint import pprint #pprint(dir(reader)) data = reader.GetOutput() print data I don't see any arrays, it implicitly says "Number Of Arrays: 0" How do I access above-mentioned data? Sincerely, George

2 2

Re: [firedrake] how to set up random initial conditions?
by Colin Cotter 26 Apr '16

26 Apr '16

Hi Francis, Yes, if you have a Firedrake function f, then f.dat.shape gives you the shape of f.dat.data, the numpy array of basis coefficients. all the best --cjc On 24 April 2016 at 00:23, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello Ali, > > Thanks a lot. That is very helpful. It is very nice how you can access > numpy functions through this syntax. Very clean. > > I tired what you suggested and I believe that the example you suggested > shifted the whole solution by a random constant. If I wanted to have a > different random variable at every position, then i presume I need > > C=np.random.normal(0,1,len(something)) > > where something is the size of the mesh. In a previous email David > suggested using > > *f.dat.shape > > I presume that f.dat is based on the mesh but I don't know how to extract > that information from the mesh. > > Again, sorry for not being able to figure this out on my own but your help > is greatly appreciated. > > Francis > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] > on behalf of Gregory, Alastair C A [a.gregory14(a)imperial.ac.uk] > *Sent:* Saturday, April 23, 2016 6:42 PM > > *To:* firedrake > *Subject:* Re: [firedrake] how to set up random initial conditions? > > Hi Francis, > > > Thanks for the clarification. The f.dat.data just comes from accessing the > basis coefficients in a different way to interpolating an expression. > > If you want to input the randomness via the expression instance then you > can manipulate the code in the following way depending on what you want to. > For example if you wanted to shift the initial condition in the X-axis by > the normally distributed variable 'C' you do the following : > > q0=Function(Vdg).interpolate(Expression("0.1*(sin(x[0])+C)*sin(x[1]) > ",C=np.random.normal(0,1,1))) > > I hope this is what you're looking for. > > Many Thanks, > Ali > > Sent from my iPhone > > On 23 Apr 2016, at 23:23, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > > Hello, > > Thanks David and Alastair for the advice. > > I am sorry that I should have been more precise in terms of my confusion. > > If you wanted to see the same code that I put together (with lots of help > from Colin I should add) you can go to the following, > > > https://github.com/francispoulin/firedrake/blob/qg_demo/demos/quasigeostrop… > > Below are the lines where we define the function space and then define the > initial conditions. > > Vdg = FunctionSpace(mesh,"DG",1) # DG elements for Potential Vorticity (PV) > # Intial Conditions for PV > q0 = Function(Vdg).interpolate(Expression("0.1*sin(x[0])*sin(x[1])")) > > > I am really confused as to where the f.dat comes in. Sorry. > > I thought that if I wanted to define the initial conditions I should do it > in terms of an expression but I don't imagine that even if I import random > from numpy that random will be understood in expression. > > I guess my more specific question is how do I use random in an expression? > > Cheers, Francis > > > ------------------ > Francis Poulin > Associate Professor > Department of Applied Mathematics > University of Waterloo > > email: fpoulin(a)uwaterloo.ca > Web: https://uwaterloo.ca/poulin-research-group/ > Telephone: +1 519 888 4567 x32637 > > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk [firedrake-bounces(a)imperial.ac.uk] > on behalf of Gregory, Alastair C A [a.gregory14(a)imperial.ac.uk] > *Sent:* Saturday, April 23, 2016 6:53 AM > *To:* firedrake > *Subject:* Re: [firedrake] how to set up random initial conditions? > > Hi Francis and David, > > From experience in using Firedrake and Uncertainty Quantification, I have > found that using spatially correlated random initial condition for a field > is better for some cases. > > You can do this by solving a simple Helmholtz equation once with the > forcing f as a random field, > > Np.random.normal(0,1,len(F.dat.data)) > > > Many Thanks, > Alastair Gregory > > Sent from my iPhone > > On 23 Apr 2016, at 11:32, David Ham <David.Ham(a)imperial.ac.uk> wrote: > > 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 >> >> >> _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > _______________________________________________ > 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 4

parallelization
by Francis Poulin 25 Apr '16

25 Apr '16

Hello, So far I have only been running firedrake runs using python job.py and I thought it was in serial. I just tried a run with a higher resolution on my dell that has 2x8 cores and it seems to be using nearly 32,000% of my cpu. I am very happy to know that it is parallelizing automatically, it makes my life a lot easier but I have two questions. 1) What kind of parallelization is it using? From the high percentage it seems like it's threaded but not sure how to check. 2) Is there any advantage to running mpirun -np 16 python job.py versus doing the more simple execution above? If not I will continue doing the simple approach. Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin(a)uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637

1 0

Re: [firedrake] Symmetry axis boundary condition
by David Ham 25 Apr '16

25 Apr '16

Hi Fryderyk, Pasting Python code into emails is always unsafe. In this case, the indentation has been mucked up, which changes the meaning of the code. Can you please share a git repo (eg on bitbucket or github) with the code so I can try running it. Regards, David On Sun, 24 Apr 2016 at 13:47 Fryderyk Wilczynski <scfw(a)leeds.ac.uk> wrote: > 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)) > >

2 1

use of BCs for interior nodes
by Floriane Gidel [RPG] 25 Apr '16

25 Apr '16

Dear all, In the code attached, the objective is to copy the surface values of the 3D function phi_s, to the interior nodes of the 3D function phi_3s such that phi_3s is equal to the top values of phi_2s at each z. For this, I want to solve dzz(phi) = 0, with boundary conditions : phi_3s = phi_s(top) at the bottom z=0 and top z=H0. solving the equation: phi.dx(2)*v.dx(2) = 0 (where phi is the trial function and v the test function) did not work, so instead I solve: grad(phi).grad(v) = phi_3s.dx(0)v.dx(0) + phi_3s.dx(1) v.dx(1), with phi_s as boundary condition. I have two questions concerning this code: 1) First, when I change CG1 to CG2, then the value of phi_3s is not exactly the same as the one of phi_s at the surface. Where does this difference come from? 2) Is there a way to define a boundary condition at the bottom of the domain V (at z=0), using the value of a function at the top of this domain? Basically, I would like to use the top surface value of phi_s as a boundary condition both at the bottom z=0 and at the top z=h0 of the domain. Is that possible? Thanks and best regards, Floriane

3 12

Re: [firedrake] making a vertical slice with u,v,w in one field
by Colin Cotter 25 Apr '16

25 Apr '16

Hi Lawrence, Thanks for the reply. It sounds like going with (1) would be the best plan then. I'll make an issue so I can look more closely at what needs to be done. cheers --cjc On 25 April 2016 at 05:12, Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > Hi Colin, > > On 24/04/16 00:28, Colin Cotter wrote: > > 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. > > So you would like a mesh, I'll draw quads: > > x---o > | | > | | > y---p > > Where the vertical edges and the vertices are periodic (i.e. x and o > are topologically the same and similarly y and p). > > I think you might be able to do this if you talk directly to the > dmplex API. First you make the topological cells you want, then you > build a coordinate section (a functionspace) on that DM and set up > your coordinates appropriately. We'd then have to do a little bit of > plumbing to make sure that your could tell firedrake that the > coordinate field petsc provides is already discontinuous. > > > 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. > > This would imply that different basis functions in your element pull > back differently. That is Right Out at the moment. If you would > *only* like to treat these parts of the space as (u, v, w) when > solving you could do the following. Make a 3-variable mixed space: > > W = XZ * Y * P > > In your solver options do: > > 'pc_fieldsplit_0_fields': '0, 1', > 'pc_fieldsplit_1_fields': '2', > # normal fieldsplit options for 0 and 1 where 0 now refers to the > # XZY space and 1 to the P space. > > You have to turn off nested matrices for this to work. > > > 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? > > Basis functions with different pullback is difficult I think, and > requires serious re-engineering all through the form compiler > pipeline. The solver trick I just showed should "work right now". > > Cheers, > > Lawrence > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0