May 2016 - firedrake - imperial.ac.uk

Detecting NaNs
by Gregory, Alastair C A 31 May '16

31 May '16

Hi Firedrakers, I have two functions, v(x,y) and h(x,y). I find a quantity u(x,y)=v(x,y)/h(x,y) to solve a PDE. Now, I'm struggling with using u(x,y) when h(x,y) is equal to 0 at some cells. I get DIVERGED_NAN errors as I'm trying to divide by 0. Now what I have done is to build a kernel that allows me to set the cells that will be dividing by 0, to 0, and otherwise proceed as normal on the cells that won't be dividing by 0. This arises since v is momentum and h is depth, and thus u is velocity. With 0 depth, I would like 0 velocity. However, sometimes I divide by a restriction to a function and thus can't use a kernel. I was wondering if there was a clean way (potentially UFL way) to set cells of functions from NANS to 0, or anything similar that would help me? Many Thanks, Alastair Gregory ------------------------------------------------------------------------------------------ Research Assistant (Maths Foresees Feasibility Study) Imperial College London Room 759 Huxley Building, South Kensington (Tel: 07794 243913) | (Email: a.gregory14(a)imperial.ac.uk) -------------------------------------------------------------------------------------------

2 1

Two coupled problems
by Daniel Ruprecht 31 May '16

31 May '16

Dear Firedrakers, I am trying to solve two coupled problems on different domains (essentially a small box next to a larger box), where the solution at the boundary from one problem determines the Neumann boundary condition for the second problem at the interface. I am not sure how to efficiently use the solution of the first problem to prescribe the boundary condition for the second. My idea was to derive from the Expression class and overwrite the eval function including a couple of if statements to make sure the coordinates are in the region where both meshes actually match. Is this the correct way to do it (it seems quite slow, probably because the derived expression can't be compiled in C) or is there a cleverer way to do this? 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

2 6

Re: [firedrake] Caught signal number 15 Terminate
by David Ham 30 May '16

30 May '16

Alternatively, might this be the PETSc memory leak? On Mon, 30 May 2016 at 15:18 Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > > > On 30 May 2016, at 13:20, Fryderyk Wilczynski <scfw(a)leeds.ac.uk> wrote: > > > > Dear Firedrakes, > > > > My computations get killed with the following exception: > > > > > > COFFEE finished in 0.0201869 seconds (flops: 92 -> 92) > > COFFEE finished in 0.0305271 seconds (flops: 92 -> 92) > > COFFEE finished in 0.0411251 seconds (flops: 92 -> 92) > > [1] pyop2:INFO Solving nonlinear variational problem... > > [2] pyop2:INFO Solving nonlinear variational problem... > > ... > > > [1]PETSC ERROR: to get more information on the crash. > > > -------------------------------------------------------------------------- > > mpirun noticed that process rank 3 with PID 6682 on node comp-pc3234 > exited on signal 9 (Killed). > > Probably this happens because of a segfault somewhere in the code. Are > you using Expressions in your code that maybe are not safe (overwriting > memory somewhere). Does this happen in serial? It is hard to know more > without knowing exactly what you are doing. > > Lawrence > >

2 1

Re: [firedrake] saving vtu files in binary
by Colin Cotter 30 May '16

30 May '16

Don't forget that you are using P1dg data so this is 3 DOFs per cell. cheers --cjc On 30 May 2016 at 15:09, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello Lawrence, > > In a case with 100x100x100, each field seems to be taking up about 200 > MB. I talked to someone who is saving things in binary and they thought > this was probably too big by a factor of 50 or so. I'm not sure how much > space it should take but if it's possible to be more compact that would be > nice. I can compress things afterwards but in my experience binary files > usually don't compress that much, but I could be wrong. > > 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 Lawrence Mitchell [lawrence.mitchell(a)imperial.ac.uk] > Sent: Monday, May 30, 2016 9:54 AM > To: firedrake(a)imperial.ac.uk > Subject: Re: [firedrake] saving vtu files in binary > > Dear Francis, > > > On 30 May 2016, at 14:43, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > > > > Hello, > > > > I have noticed that the vtu files that I'm creating in firedrake tend to > be big. I presume this is because they are being saved in ascii format. > Is it easy to save these files in binary to be a bit more space efficient? > > Actually, only the metadata is saved in ascii format. The actual bulk > data is saved in binary (although, to be fair, it's uncompressed). > Unfortunately due to limitations in the vtu file format, the mesh > description must be saved in every file (even when it doesn't change) which > makes for some size increases. Are the vtus you get much larger than > you're expecting? Or is is just that you'd like them to be smaller? > > Cheers, > > Lawrence > > _______________________________________________ > 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

1 0

saving vtu files in binary
by Francis Poulin 30 May '16

30 May '16

Hello, I have noticed that the vtu files that I'm creating in firedrake tend to be big. I presume this is because they are being saved in ascii format. Is it easy to save these files in binary to be a bit more space efficient? 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

2 3

Re: [firedrake] saving vtu files in binary
by Andrew McRae 30 May '16

30 May '16

What function space is the field that you're saving? If it's not fully continuous then we output the data as fully discontinuous, which might have some redundancy. [I'll remark that 100*100*100*(8 bytes/double) is already ~8MB, so that's only a factor of 25. Since the mesh data is also there, multiply that 8MB by (1 + 3).] On 30 May 2016 at 15:09, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > Hello Lawrence, > > In a case with 100x100x100, each field seems to be taking up about 200 > MB. I talked to someone who is saving things in binary and they thought > this was probably too big by a factor of 50 or so. I'm not sure how much > space it should take but if it's possible to be more compact that would be > nice. I can compress things afterwards but in my experience binary files > usually don't compress that much, but I could be wrong. > > 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 Lawrence Mitchell [lawrence.mitchell(a)imperial.ac.uk] > Sent: Monday, May 30, 2016 9:54 AM > To: firedrake(a)imperial.ac.uk > Subject: Re: [firedrake] saving vtu files in binary > > Dear Francis, > > > On 30 May 2016, at 14:43, Francis Poulin <fpoulin(a)uwaterloo.ca> wrote: > > > > Hello, > > > > I have noticed that the vtu files that I'm creating in firedrake tend to > be big. I presume this is because they are being saved in ascii format. > Is it easy to save these files in binary to be a bit more space efficient? > > Actually, only the metadata is saved in ascii format. The actual bulk > data is saved in binary (although, to be fair, it's uncompressed). > Unfortunately due to limitations in the vtu file format, the mesh > description must be saved in every file (even when it doesn't change) which > makes for some size increases. Are the vtus you get much larger than > you're expecting? Or is is just that you'd like them to be smaller? > > Cheers, > > Lawrence > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >

1 0

Caught signal number 15 Terminate
by Fryderyk Wilczynski 30 May '16

30 May '16

Dear Firedrakes, My computations get killed with the following exception: COFFEE finished in 0.0201869 seconds (flops: 92 -> 92) COFFEE finished in 0.0305271 seconds (flops: 92 -> 92) COFFEE finished in 0.0411251 seconds (flops: 92 -> 92) [1] pyop2:INFO Solving nonlinear variational problem... [2] pyop2:INFO Solving nonlinear variational problem... [0] pyop2:INFO Solving nonlinear variational problem... [3] pyop2:INFO Solving nonlinear variational problem... [0] pyop2:INFO Compiling wrapper... [0] pyop2:INFO Compiling wrapper...done [2]PETSC ERROR: ------------------------------------------------------------------------ [2]PETSC ERROR: Caught signal number 15 Terminate: Some process (or the batch system) has told this process to end [2]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger [2]PETSC ERROR: or see http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind [2]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors [2]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run [2]PETSC ERROR: to get more information on the crash. [0]PETSC ERROR: ------------------------------------------------------------------------ [0]PETSC ERROR: Caught signal number 15 Terminate: Some process (or the batch system) has told this process to end [0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger [0]PETSC ERROR: or see http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind [0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors [0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run [0]PETSC ERROR: to get more information on the crash. [1]PETSC ERROR: ------------------------------------------------------------------------ [1]PETSC ERROR: Caught signal number 15 Terminate: Some process (or the batch system) has told this process to end [1]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger [1]PETSC ERROR: or see http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind [1]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors [1]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run [1]PETSC ERROR: to get more information on the crash. -------------------------------------------------------------------------- mpirun noticed that process rank 3 with PID 6682 on node comp-pc3234 exited on signal 9 (Killed). -------------------------------------------------------------------------- why does this happen/can it be avoided? Best, Fryderyk

2 1

Re: [firedrake] Laplace equation on subdomain
by David Ham 30 May '16

30 May '16

Hi Tomasz, There are two things here. The first one is cell subdomains. There is a pull request in at PETSc right now which will make this work to some extent. We'll let you know once the PETSc people accept that pull. The other issue is what you can pass to DirichletBC. That is always a boundary marker, not a cell subdomain marker - so the precise syntax you used there won't quite work. The functionality will hopefully be there soon though. On Tue, 24 May 2016 at 09:44 Tomasz Salwa [RPG] <mmtjs(a)leeds.ac.uk> wrote: > Dear all, > > > Can I invoke in Firedrake a subdomain id in a mesh generated by gmsh? > > It would be useful for me in constructing indicator functions. > > The line: > > ... > > bc = DirichletBC( V, 1, 101 ) > > ... > > causes an error: > > LookupError: 101 is not a valid marker > > although 101 is defined in gmsh as a Physical Volume. I had no problems > with Physical Surface in the same context. > > > 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:* 27 April 2016 17:32:56 > > *To:* firedrake > *Subject:* Re: [firedrake] Laplace equation on subdomain > > > That did the trick, I have a coupled solver working now, thank you! > > 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:* 26 April 2016 22:29:44 > *To:* firedrake > *Subject:* Re: [firedrake] Laplace equation on subdomain > Well the following reproduces the reported error for me: > > from firedrake import * > m = UnitSquareMesh(2, 2) > fs = FunctionSpace(m, "CG", 1) > v = TestFunction(fs) > assemble(v*dS) > > > On Tue, 26 Apr 2016 at 21:27 Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: > >> I didn't think this was true, unless something changed "recently". >> >> I can't test right now, but I'm quite confident I used to use continuous >> functions in interior facet integrals on simplex meshes without needing a >> restriction. >> >> However, I recall being forced to use restrictions on extruded meshes; I >> never got to the bottom of why. Perhaps Miklós can comment (or state that >> my recollection is nonsense)? >> >> On 26 April 2016 at 21:12, David Ham <David.Ham(a)imperial.ac.uk> wrote: >> >>> 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 >>>>>>> >>>>>>

3 3

Product vector matrix
by Floriane Gidel [RPG] 27 May '16

27 May '16

Dear all, I would like to solve equations in which one of the unknowns is a vector that I need to multiply to some matrices. However, these matrices are not computed with the assemble function but with a python function and are thus arrays. To put this into context, I am solving 3D equations with a separation of variable between the (x,y) dependency and the z dependency. Then, the z-discretization is made "by hand" and then substituted into the horizontal weak formulations which are solved with Firedrake. Therefore, the discrete vertical matrices are computed separately, but I would like to include them into my weak formulations. The issue is that it does not seem to work when multiplying a function with an array in the variational problem. How could I do that ? The WF and a description of the objective is in the attached PDF. Thanks in advance, Floriane

1 0

Re: [firedrake] Applying DirichletBCs to arbitrary \Gamma_{inflow} regions
by Colin Cotter 26 May '16

26 May '16

This all seems a bit odd. Inflow boundary conditions for the advection equation should be handled weakly, then you can make an upwind flag to do this as normal. cheers --cjc On 26 May 2016 at 09:19, Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > Hi Justin, > > > On 26 May 2016, at 08:43, Justin Chang <jychang48(a)gmail.com> wrote: > > > > Hi all, > > > > I am attempting to simulate a chaotic advection equation where the > velocity vector field is based on the Arnold–Beltrami–Childress (ABC) flow > equation: > > > > v_x = A*sin(z) + C*cos(y) > > v_x = B*sin(x) + A*cos(z) > > v_x = C*sin(y) + B*cos(x) > > > > Right now I have the following: > > > > from firedrake import * > > mesh = UnitCubeMesh(20,20,20) > > A = 0.08 > > B = 0.1 > > C = 0.15 > > V = VectorFunctionSpace(mesh, "CG", 1) > > velocity = Function(V) > > velocity.interpolate(Expression(( > > "A*sin(2*pi*x[2])+C*cos(2*pi*x[1])", > > "B*sin(2*pi*x[0])+A*cos(2*pi*x[2])", > > "C*sin(8*pi*x[1])+B*cos(8*pi*x[0])"), > > A=A,B=B,C=C)) > > FWIW, I recommend for interpolating expressions to utilise the new UFL > expression interpolation. > http://firedrakeproject.org/interpolation.html#ufl-expressions This gets > you much better type checking and more useful error messages in the case of > typos. > > > Where A,B,C are user-defined constants (will add such constants inside > all sin and cos functions later). I want to have DirichletBC's only applied > to faces and regions where the velocity points into the cube (i.e., the > "upstream" location). > > This is doable, you need to subclass DirichletBC and provide your own > property that computes the *nodes* at which the boundary condition should > apply. > > In the following, I assume you want to apply a bc to the same space that > the velocity lives in. I would do the following. > > import operator > import numpy as np > > > class UpstreamBC(DirichletBC): > def __init__(self, V, value, sub_domain, method="topological", > velocity=None): > super(UpstreamBC, self).__init__(V, value, sub_domain, > method=method) > self._velocity = velocity > self._all_nodes = super(UpstreamBC, self).nodes > > @property > def nodes(self): > face_velocity = > self._velocity.dat.data_ro_with_halos[self._all_nodes] > component, comparator = {1: (0, operator.gt), # x == 0 > 2: (0, operator.lt), # x == 1 > 3: (1, operator.gt), # y == 0 > 4: (1, operator.lt), # y == 1 > 5: (2, operator.gt), # z == 0 > 6: (2, operator.lt), # z == 1 > }[self.sub_domain] > return self._all_nodes[comparator(face_velocity[:, component], > 0.0)] > > @property > def node_set(self): > return op2.Subset(self._function_space.node_set, self.nodes) > > Untested, but you can see the idea. I provide, in addition the the normal > arguments an additional velocity field. In the nodes property I use the > value of the velocity field to select the appropriate subset of the nodes. > > Note that to make this selection simple, I require that the bc object be > made with a single subdomain value (that way I can simply determine which > component of the velocity field I should be checking). So you need one BC > object for each side of the cube. This is possibly slightly inefficient, > so we could improve it, but let's try and get everything working first. > > > The region locations and sizes will vary depending on the constants > chosen. Yes I could create separate GMSH files and manually determine where > these "inflow" regions will be for every set of A/B/C/etc constants, but I > was wondering if there's a way to automate this the Firedrake way, > particularly if these constants will vary as a function of time. > > The setup I have above is already time-dependent, so you don't need to do > anything else. > > Hope this works! > > Lawrence > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

2 1