August 2015 - firedrake - imperial.ac.uk

Multi-grid preconditioning
by Justin Chang 16 Nov '15

16 Nov '15

Hi everyone, Is there support for Multi-grid preconditioning? I think I saw somewhere where this was not supported yet. Because applying -pc_type gamg for my mixed forms returns me errors. Thanks, Justin

2 3

Installing Firedrake on an HPC machine
by Justin Chang 14 Nov '15

14 Nov '15

Hi everyone, I have installed firedrake on my university's HPC machine, and whenever i attempt to run any Firedrake program, I get this error: -------------------------------------------------------------------------- An MPI process has executed an operation involving a call to the "fork()" system call to create a child process. Open MPI is currently operating in a condition that could result in memory corruption or other system errors; your MPI job may hang, crash, or produce silent data corruption. The use of fork() (or system() or other calls that create child processes) is strongly discouraged. The process that invoked fork was: Local host: compute-0-0 (PID 28214) MPI_COMM_WORLD rank: 0 If you are *absolutely sure* that your application will successfully and correctly survive a call to fork(), you may disable this warning by setting the mpi_warn_on_fork MCA parameter to 0. -------------------------------------------------------------------------- Why is this happening? Program seems fine when I only use a few processes but when I use two or more compute nodes, my program hangs. I'll do my best to explain how I built firedrake: 1) Our system has openmpi-1.8.3 and gcc-4.9.2 available. so those are the only modules I have loaded, everything else was downloaded. 2) I built Python 2.7.10 from source and installed it in ${HOME}/.local 3) pip was obtained via python get-pip.py 4) PETSc was installed --with-cc=mpmicc, --with-cxx=mpicxx, and has all the necessary external packages 5) petsc4py was installed via git clone ... and set up with python setup.py build/install 6) PyOP2 was built with python setup.py build_ext --inplace 7) Firedrake was built via make and everything else was installed via 'pip install <package>'. I had to manually install swig and pcre. Any ideas what this may be? And/or how I can circumvent this issue? Thanks, Justin

2 20

Re: [firedrake] Solve a Variational problem in a part of the domain
by Colin Cotter 30 Sep '15

30 Sep '15

Well, can you sketch out how you would solve this system, then we can see how to build it. all the best --cjc On 17 August 2015 at 13:33, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> wrote: > Yes, obviously I wouldn't be able to write > > assemble(u*dx) > > for the trial function u, but I just wanted to show you what I have. > > Is it too difficult to design such a non-conventional solver? > > Thanks, > > Anna. > > > > On 17/08/15 13:04, Colin Cotter wrote: > > Hi Anna, > Not as written, no. I think that we would need to design a bespoke > solver for that. > > all the best > --cjc > > On 17 August 2015 at 12:18, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> > wrote: > >> Hi all, >> >> I have to solve a problem with the bilinear form similar to the following: >> >> a = ( v*u + dt*inner(grad(u),grad(v)) )*dx + ( v*dt*assemble(u*dx) )*dx >> >> The last term is essentially the product of the integral of u and the >> integral of test function v. >> Is it even possible to solve this in Firedrake? >> >> Best, >> >> Anna. >> >> >> On 14/08/15 15:43, Lawrence Mitchell wrote: >> >> On 14 Aug 2015, at 11:55, Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> <a.kalogirou(a)leeds.ac.uk> wrote: >> >> Hi, >> >> I was planning to define the heavyside function as a DG0 function. >> >> Any ideas about solving equation (1b), which contains a time-update for both a function and a scalar? >> >> I haven't been fully following along, but it looks like the formulation in (1f) introduces, effectively, a global coupling. Is that right? Assuming you didn't have any implementation constraints, how would you go about solving this? >> >> Lawrence >> >> >> >> _______________________________________________ >> firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >> -- >> >> Dr Anna Kalogirou >> Research Fellow >> School of Mathematics >> University of Leeds >> >> http://www1.maths.leeds.ac.uk/~matak/ >> >> > > > -- > http://www.imperial.ac.uk/people/colin.cotter > > www.cambridge.org/9781107663916 > > > > > _______________________________________________ > firedrake mailing listfiredrake@imperial.ac.ukhttps://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > -- > > Dr Anna Kalogirou > Research Fellow > School of Mathematics > University of Leeds > > http://www1.maths.leeds.ac.uk/~matak/ > > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

3 3

Few questions
by Justin Chang 01 Sep '15

01 Sep '15

Hi all, 1) How do I define a volumetric source/sink for a specific region? For example, if I am solving a diffusion equation with homogeneous BCs and have a volumetric source f(x,y) such that if (x[0] >= 0.25 && x[0] <= 0.30 && x[1] >= 0.25 && x[1] <= 0.30): f == 1 else: f == 0 2) I want to play around with PETSc's Variational Inequality <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/SNES/SNESVINEWTONSS…> solver. In my own PETSc code all I had to do was provide a Jacobian matrix, a residual vector, pass in the options -snes_type vinewtonssls (or it might be vinewotnrsls), and pass in SNESVINEWTONSSLS <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/SNES/SNESVINEWTONSS…> . >From my last discussion with some PETSc folks it seems this VI feature is relatively now, so I was wondering if that last item can be invoked through the version of petsc4py you guys are using? 3) Is it possible to employ the SUPG formulation for advection-diffusion problems in firedrake? Thanks, Justin

2 3

Re: [firedrake] Few questions
by Andrew McRae 31 Aug '15

31 Aug '15

1) Does this region line up with cells of your mesh? If so, my normal way is to interpolate that expression into a P0 function: P0 = Function(mesh, "DG", 0) f = Function(P0) f.interpolate(Expression("(x[0] >= 0.25 && x[0] <= 0.30 && x[1] >= 0.25 && x[1] <= 0.30) ? 1.0 : 0.0")) and then use f in calculations. 2) Pass, someone else can answer this 3) "Yes", but can you be more specific? SUPG is just modifying the test function, so my code looked something like ... u = ... alpha = ... phi = TestFunction(V) phiSUPG = phi + alpha*dot(u, grad(phi)) a = phiSUPG*....*dx + .... On 31 August 2015 at 22:31, Justin Chang <jychang48(a)gmail.com> wrote: > Hi all, > > 1) How do I define a volumetric source/sink for a specific region? For > example, if I am solving a diffusion equation with homogeneous BCs and have > a volumetric source f(x,y) such that > > if (x[0] >= 0.25 && x[0] <= 0.30 && x[1] >= 0.25 && x[1] <= 0.30): > f == 1 > else: > f == 0 > > 2) I want to play around with PETSc's Variational Inequality > <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/SNES/SNESVINEWTONSS…> > solver. In my own PETSc code all I had to do was provide a Jacobian matrix, > a residual vector, pass in the options -snes_type vinewtonssls (or it might > be vinewotnrsls), and pass in SNESVINEWTONSSLS > <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/SNES/SNESVINEWTONSS…> . > From my last discussion with some PETSc folks it seems this VI feature is > relatively now, so I was wondering if that last item can be invoked through > the version of petsc4py you guys are using? > > 3) Is it possible to employ the SUPG formulation for advection-diffusion > problems in firedrake? > > Thanks, > Justin >

1 0

Re: [firedrake] The wonderful world of firedrake-install
by David Ham 30 Aug '15

30 Aug '15

Hi All, I have a branch going through testing now which will handle PETSC_CONFIGURE_OPTIONS in a smarter way. Basically we'll just make sure that the things which are required are there, and we'll also honour anything the user has already set. I've also updated docs to say that that's what happens. I am also open to the suggestion that we could add more configuration options to the default set. Would anyone like to suggest what they should be? Cheers, David On Fri, 28 Aug 2015 at 10:25 Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > > > On 27 Aug 2015, at 23:57, Justin Chang <jychang48(a)gmail.com> wrote: > > > > David, > > > > A few comments: > > > > 1) If I run that script (w/ --developer) and have the latest installed > > petsc-dev and established $PETSC_DIR, the installation of petsc4py > > will produce an error similar to the one reported in this thread: > > http://lists.mcs.anl.gov/pipermail/petsc-users/2015-July/026404.html > > to which I had to either: a) use the patch Satish provided or b) clone > > from https://bitbucket.org/petsc/petsc4py and not mapdes/petsc4py. If > > I unset $PETSC_DIR then the installation works fine. Not a major > > problem or anything but it seems to me that the installation script is > > using a somewhat outdated petsc4py so I was wondering if this was > > intentional or not. > > > The installation uses a matching set of PETSc/petsc4py, that are > "up-to-date" and known good. But unfortunately it appears that if > PETSC_DIR is set to a more recent PETSc, then the downloaded petsc4py might > not match. The script provides a warning in this case, but it was maybe > missed in all the other output. Perhaps if PETSC_DIR is set the script > should exit with an error (overridable if you definitely want to use the > external PETSc). > > > 2) Some of my problems require some PETSc external packages like > > hypre, but I see that the default configure options only include > > chaco, ctegen, and triangle. I had to manually add the > > --download-hypre to the firedrake-install to get what I wanted, but do > > you think y'all can also include other important packages like mumps, > > superlu_dist, hdf5, metis, parmetis, etc? > > You can add options to the petsc configure call by setting the environment > variable PETSC_CONFIGURE_OPTIONS to appropriate values. Make sure you > include the necessary packages (chaco, ctetgen and triangle). Perhaps > instead one should be able specify additional PETSc configure options with > something like: > > firedrake-install ... --petsc-options="..." > > > Cheers, > > Lawrence >

3 2

The wonderful world of firedrake-install
by David Ham 28 Aug '15

28 Aug '15

Dear Firedrakers, I am happy to be able to announce the release of our new install and upgrade system for Firedrake. The new system makes installing and upgrading Firedrake much easier than the previous instruction set. Particular features: * single command install on MacOS X, Ubuntu (and sometimes other Linux). * support for per user or global installations. * support for installation into a Python virtualenv to ensure no interference between Firedrake and other software (in particular this makes installing Firedrake alongside Dolfin safer and easier). * Support for "developer mode" installs where core Firedrake components run from the source directories. The core Firedrake developers will be using this install, and it's the install system used by our automated build testing, so you can be assured that it will always be well tested and emerging bugs will be fixed fast. Installation instructions are available now from http://www.firedrakeproject.org/download.html I look forward to your feedback on the new system. Cheers, David

3 2

Re: [firedrake] Specify boundary conditions at a corner
by David Ham 28 Aug '15

28 Aug '15

Hi Justin, Thanks for that, you uncovered a bug in some of how our internal imports work. I have just merged a fix into Firedrake master so hopefully this should now work. Cheers, David On Wed, 26 Aug 2015 at 10:29 Justin Chang <jychang48(a)gmail.com> wrote: > David, > > I am still getting the same error. Attached the code, can you have a > look at it and see what's going on? Note that if you comment out the > point dirichlet stuff and uncomment the nullspace related operations, > it outputs the correct solution. > > Thanks, > Justin > > On Tue, Aug 25, 2015 at 3:03 AM, David Ham <David.Ham(a)imperial.ac.uk> > wrote: > > Hi Justin, > > > > The imports in bcs.py are relative imports (which is slightly bad style > > but...) In order to access correct module from outside you need to use: > > > > from firedrake import utils > > > > Regards, > > > > David > > > > On Tue, 25 Aug 2015 at 00:48 Justin Chang <jychang48(a)gmail.com> wrote: > >> > >> David, > >> > >> I get the following error: > >> > >> Traceback (most recent call last): > >> File "Test_LSAug.py", line 41, in <module> > >> class PointDirichletBC(DirichletBC): > >> File "Test_LSAug.py", line 42, in PointDirichletBC > >> @utils.cached_property > >> AttributeError: 'module' object has no attribute 'cached_property' > >> > >> When I comment out that line, i get a long error. When I try to import > >> utils and all the other modules that bcs.py imports, i get errors > >> saying those modules don't exist. Know what's going on? > >> > >> Thanks, > >> Justin > >> > >> On Mon, Aug 24, 2015 at 3:14 AM, David Ham <David.Ham(a)imperial.ac.uk> > >> wrote: > >> > Hi Justin, > >> > > >> > Having subclassed DirichletBC, what you now have is a new type of > >> > Dirichlet > >> > "boundary" object, so you just use it as a Dirichlet condition: > >> > > >> > bc3 = PointDirichletBC(W.sub(1), 0, 0) > >> > > >> > bc_all = [bc1, bc2, bc3] > >> > > >> > However Andrew is right, removing the nullspace is almost certainly a > >> > better > >> > solution in this case. > >> > > >> > Cheers, > >> > > >> > David > >> > > >> > On Sun, 23 Aug 2015 at 05:23 Justin Chang <jychang48(a)gmail.com> > wrote: > >> >> > >> >> Andrew, > >> >> > >> >> That actually did the trick, thank you very much. > >> >> > >> >> But I still would like to know the answer to my original question. Or > >> >> perhaps, instead of a single corner, I would like to set (x=0,y<0.1) > >> >> and (x<0.1,y=0) to a specific value while everywhere else is set to a > >> >> different value. > >> >> > >> >> Thanks, > >> >> Justin > >> >> > >> >> On Sat, Aug 22, 2015 at 11:58 AM, Andrew McRae > >> >> <a.mcrae12(a)imperial.ac.uk> > >> >> wrote: > >> >> > This doesn't answer your main question, but setting an appropriate > >> >> > nullspace > >> >> > might be more appropriate than pinning a single value; see > >> >> > http://www.firedrakeproject.org/solving-interface.html > >> >> > > >> >> > [If I'm wrong, I'll let the usual suspects correct my > misinformation] > >> >> > > >> >> > On 22 August 2015 at 18:47, Justin Chang <jychang48(a)gmail.com> > wrote: > >> >> >> > >> >> >> David, > >> >> >> > >> >> >> How exactly do I use that class in my code? Say I have the > following > >> >> >> function spaces/discretization based on the least-squares finite > >> >> >> element method: > >> >> >> > >> >> >> # Mesh > >> >> >> mesh = UnitSquareMesh(seed,seed) > >> >> >> V = VectorFunctionSpace(mesh,"CG",1) > >> >> >> Q = FunctionSpace(mesh,"CG",1) > >> >> >> W = V*Q > >> >> >> v,p = TrialFunctions(W) > >> >> >> w,q = TestFunctions(W) > >> >> >> > >> >> >> # Weak form > >> >> >> g = Function(V) > >> >> >> > >> >> >> > >> >> >> > >> >> >> > g.interpolate(Expression(("cos(pi*x[0])*sin(pi*x[1])+2*pi*cos(2*pi*x[0])*sin(2*pi*x[1])","-sin(pi*x[0])*cos(pi*x[1])+2*pi*sin(2*pi*x[0])*cos(2*pi*x[1])"))) > >> >> >> a = dot(v+grad(p),w+grad(q))*dx + div(v)*div(w)*dx > >> >> >> L = dot(w+grad(q),g)*dx > >> >> >> > >> >> >> # Boundary conditions > >> >> >> bc1 = DirichletBC(W.sub(0).sub(0), > >> >> >> Expression("cos(pi*x[0])*sin(pi*x[1])"), (1,2)) > >> >> >> bc2 = DirichletBC(W.sub(0).sub(1), > >> >> >> Expression("-sin(pi*x[0])*cos(pi*x[1])"), (3,4)) > >> >> >> bc_all = [bc1,bc2] > >> >> >> > >> >> >> # Solve > >> >> >> cg_parameters = { > >> >> >> 'ksp_type': 'cg', > >> >> >> 'pc_type': 'bjacobi' > >> >> >> } > >> >> >> solution = Function(W) > >> >> >> A = assemble(a,bcs=bc_all) > >> >> >> b = assemble(L,bcs=bc_all) > >> >> >> solver = > >> >> >> > LinearSolver(A,solver_parameters=cg_parameters,options_prefix="cg_") > >> >> >> solver.solve(solution,b) > >> >> >> > >> >> >> If I run the above code I get an error saying 'LinearSolver failed > >> >> >> to > >> >> >> converge after %d iterations with reason: %s', 196, > >> >> >> 'DIVERGED_INDEFINITE_MAT'. Which I am guessing has to do with the > >> >> >> lack > >> >> >> of a boundary condition for the Q space, thus I want to ensure a > >> >> >> unique solution by prescribing the bottom left constraint to a > zero > >> >> >> value. > >> >> >> > >> >> >> Thanks, > >> >> >> Justin > >> >> >> > >> >> >> > >> >> >> On Fri, Aug 21, 2015 at 4:52 AM, David Ham > >> >> >> <David.Ham(a)imperial.ac.uk> > >> >> >> wrote: > >> >> >> > Hi Justin, > >> >> >> > > >> >> >> > The nice way of doing this would require better subdomain > support > >> >> >> > than > >> >> >> > we > >> >> >> > have right now. However there is a slightly hacky way of doing > it > >> >> >> > which > >> >> >> > I > >> >> >> > think will cover your case nicely. > >> >> >> > > >> >> >> > If you take a look at the DirichletBC class (in bcs.py), you'll > >> >> >> > notice > >> >> >> > that > >> >> >> > the set of nodes at which the BC should be applied is calculated > >> >> >> > in > >> >> >> > DirichletBC.nodes . So you could simply subclass DirichletBC and > >> >> >> > replace > >> >> >> > nodes with a function which returns the index of the zero node. > >> >> >> > For > >> >> >> > example > >> >> >> > (and I confess this is a sketch code which I haven't tried to > >> >> >> > run): > >> >> >> > > >> >> >> > class PointDirichletBC(DirichletBC): > >> >> >> > @utils.cached_property > >> >> >> > def nodes(self): > >> >> >> > # Find the array of coordinate values. > >> >> >> > x = self.function_space().mesh().coordinates.dat.data_ro > >> >> >> > # Find the location of the zero rows in that > >> >> >> > return np.where(~x.any(axis=1))[0] > >> >> >> > > >> >> >> > Does that work for you? > >> >> >> > > >> >> >> > Cheers, > >> >> >> > > >> >> >> > David > >> >> >> > > >> >> >> > On Fri, 21 Aug 2015 at 03:32 Justin Chang <jychang48(a)gmail.com> > >> >> >> > wrote: > >> >> >> >> > >> >> >> >> Hi all, > >> >> >> >> > >> >> >> >> If I create a mesh using UnitSquareMesh or UnitCubeMesh, is > there > >> >> >> >> a > >> >> >> >> way to subject a single point (as opposed to an entire > edge/face) > >> >> >> >> to > >> >> >> >> a > >> >> >> >> DirichletBC? I want to subject the the location x=0,y=0 to some > >> >> >> >> value. > >> >> >> >> > >> >> >> >> Thanks, > >> >> >> >> Justin > >> >> >> >> > >> >> >> >> _______________________________________________ > >> >> >> >> 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 > >> >> >> > > >> >> >> > >> >> >> _______________________________________________ > >> >> >> 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 > >> >> > > >> >> > >> >> _______________________________________________ > >> >> 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 > >> > > >> > >> _______________________________________________ > >> 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 > > >

3 2

Re: [firedrake] Specify boundary conditions at a corner
by David Ham 26 Aug '15

26 Aug '15

Hi Justin, The nice way of doing this would require better subdomain support than we have right now. However there is a slightly hacky way of doing it which I think will cover your case nicely. If you take a look at the DirichletBC class (in bcs.py), you'll notice that the set of nodes at which the BC should be applied is calculated in DirichletBC.nodes . So you could simply subclass DirichletBC and replace nodes with a function which returns the index of the zero node. For example (and I confess this is a sketch code which I haven't tried to run): class PointDirichletBC(DirichletBC): @utils.cached_property def nodes(self): # Find the array of coordinate values. x = self.function_space().mesh().coordinates.dat.data_ro # Find the location of the zero rows in that return np.where(~x.any(axis=1))[0] Does that work for you? Cheers, David On Fri, 21 Aug 2015 at 03:32 Justin Chang <jychang48(a)gmail.com> wrote: > Hi all, > > If I create a mesh using UnitSquareMesh or UnitCubeMesh, is there a > way to subject a single point (as opposed to an entire edge/face) to a > DirichletBC? I want to subject the the location x=0,y=0 to some value. > > Thanks, > Justin > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >

3 2

Re: [firedrake] Specify boundary conditions at a corner
by David Ham 26 Aug '15

26 Aug '15

Hi Justin, The imports in bcs.py are relative imports (which is slightly bad style but...) In order to access correct module from outside you need to use: from firedrake import utils Regards, David On Tue, 25 Aug 2015 at 00:48 Justin Chang <jychang48(a)gmail.com> wrote: > David, > > I get the following error: > > Traceback (most recent call last): > File "Test_LSAug.py", line 41, in <module> > class PointDirichletBC(DirichletBC): > File "Test_LSAug.py", line 42, in PointDirichletBC > @utils.cached_property > AttributeError: 'module' object has no attribute 'cached_property' > > When I comment out that line, i get a long error. When I try to import > utils and all the other modules that bcs.py imports, i get errors > saying those modules don't exist. Know what's going on? > > Thanks, > Justin > > On Mon, Aug 24, 2015 at 3:14 AM, David Ham <David.Ham(a)imperial.ac.uk> > wrote: > > Hi Justin, > > > > Having subclassed DirichletBC, what you now have is a new type of > Dirichlet > > "boundary" object, so you just use it as a Dirichlet condition: > > > > bc3 = PointDirichletBC(W.sub(1), 0, 0) > > > > bc_all = [bc1, bc2, bc3] > > > > However Andrew is right, removing the nullspace is almost certainly a > better > > solution in this case. > > > > Cheers, > > > > David > > > > On Sun, 23 Aug 2015 at 05:23 Justin Chang <jychang48(a)gmail.com> wrote: > >> > >> Andrew, > >> > >> That actually did the trick, thank you very much. > >> > >> But I still would like to know the answer to my original question. Or > >> perhaps, instead of a single corner, I would like to set (x=0,y<0.1) > >> and (x<0.1,y=0) to a specific value while everywhere else is set to a > >> different value. > >> > >> Thanks, > >> Justin > >> > >> On Sat, Aug 22, 2015 at 11:58 AM, Andrew McRae < > a.mcrae12(a)imperial.ac.uk> > >> wrote: > >> > This doesn't answer your main question, but setting an appropriate > >> > nullspace > >> > might be more appropriate than pinning a single value; see > >> > http://www.firedrakeproject.org/solving-interface.html > >> > > >> > [If I'm wrong, I'll let the usual suspects correct my misinformation] > >> > > >> > On 22 August 2015 at 18:47, Justin Chang <jychang48(a)gmail.com> wrote: > >> >> > >> >> David, > >> >> > >> >> How exactly do I use that class in my code? Say I have the following > >> >> function spaces/discretization based on the least-squares finite > >> >> element method: > >> >> > >> >> # Mesh > >> >> mesh = UnitSquareMesh(seed,seed) > >> >> V = VectorFunctionSpace(mesh,"CG",1) > >> >> Q = FunctionSpace(mesh,"CG",1) > >> >> W = V*Q > >> >> v,p = TrialFunctions(W) > >> >> w,q = TestFunctions(W) > >> >> > >> >> # Weak form > >> >> g = Function(V) > >> >> > >> >> > >> >> > g.interpolate(Expression(("cos(pi*x[0])*sin(pi*x[1])+2*pi*cos(2*pi*x[0])*sin(2*pi*x[1])","-sin(pi*x[0])*cos(pi*x[1])+2*pi*sin(2*pi*x[0])*cos(2*pi*x[1])"))) > >> >> a = dot(v+grad(p),w+grad(q))*dx + div(v)*div(w)*dx > >> >> L = dot(w+grad(q),g)*dx > >> >> > >> >> # Boundary conditions > >> >> bc1 = DirichletBC(W.sub(0).sub(0), > >> >> Expression("cos(pi*x[0])*sin(pi*x[1])"), (1,2)) > >> >> bc2 = DirichletBC(W.sub(0).sub(1), > >> >> Expression("-sin(pi*x[0])*cos(pi*x[1])"), (3,4)) > >> >> bc_all = [bc1,bc2] > >> >> > >> >> # Solve > >> >> cg_parameters = { > >> >> 'ksp_type': 'cg', > >> >> 'pc_type': 'bjacobi' > >> >> } > >> >> solution = Function(W) > >> >> A = assemble(a,bcs=bc_all) > >> >> b = assemble(L,bcs=bc_all) > >> >> solver = > >> >> LinearSolver(A,solver_parameters=cg_parameters,options_prefix="cg_") > >> >> solver.solve(solution,b) > >> >> > >> >> If I run the above code I get an error saying 'LinearSolver failed to > >> >> converge after %d iterations with reason: %s', 196, > >> >> 'DIVERGED_INDEFINITE_MAT'. Which I am guessing has to do with the > lack > >> >> of a boundary condition for the Q space, thus I want to ensure a > >> >> unique solution by prescribing the bottom left constraint to a zero > >> >> value. > >> >> > >> >> Thanks, > >> >> Justin > >> >> > >> >> > >> >> On Fri, Aug 21, 2015 at 4:52 AM, David Ham <David.Ham(a)imperial.ac.uk > > > >> >> wrote: > >> >> > Hi Justin, > >> >> > > >> >> > The nice way of doing this would require better subdomain support > >> >> > than > >> >> > we > >> >> > have right now. However there is a slightly hacky way of doing it > >> >> > which > >> >> > I > >> >> > think will cover your case nicely. > >> >> > > >> >> > If you take a look at the DirichletBC class (in bcs.py), you'll > >> >> > notice > >> >> > that > >> >> > the set of nodes at which the BC should be applied is calculated in > >> >> > DirichletBC.nodes . So you could simply subclass DirichletBC and > >> >> > replace > >> >> > nodes with a function which returns the index of the zero node. For > >> >> > example > >> >> > (and I confess this is a sketch code which I haven't tried to run): > >> >> > > >> >> > class PointDirichletBC(DirichletBC): > >> >> > @utils.cached_property > >> >> > def nodes(self): > >> >> > # Find the array of coordinate values. > >> >> > x = self.function_space().mesh().coordinates.dat.data_ro > >> >> > # Find the location of the zero rows in that > >> >> > return np.where(~x.any(axis=1))[0] > >> >> > > >> >> > Does that work for you? > >> >> > > >> >> > Cheers, > >> >> > > >> >> > David > >> >> > > >> >> > On Fri, 21 Aug 2015 at 03:32 Justin Chang <jychang48(a)gmail.com> > >> >> > wrote: > >> >> >> > >> >> >> Hi all, > >> >> >> > >> >> >> If I create a mesh using UnitSquareMesh or UnitCubeMesh, is there > a > >> >> >> way to subject a single point (as opposed to an entire edge/face) > to > >> >> >> a > >> >> >> DirichletBC? I want to subject the the location x=0,y=0 to some > >> >> >> value. > >> >> >> > >> >> >> Thanks, > >> >> >> Justin > >> >> >> > >> >> >> _______________________________________________ > >> >> >> 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 > >> >> > > >> >> > >> >> _______________________________________________ > >> >> 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 > >> > > >> > >> _______________________________________________ > >> 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 > > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >

2 1