Re: [firedrake] Specify boundary conditions at a corner
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@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@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@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@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@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@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@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@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@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@imperial.ac.uk > >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > > > > _______________________________________________ > > firedrake mailing list > > firedrake@imperial.ac.uk > > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > > _______________________________________________ > firedrake mailing list > firedrake@imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
David, Okay yes this issue is resolve now, thank you very much. Lawrence, So if I use a GMSH file, the third argument in DirichletBC(...) would correspond to whatever the Physical IDS are? Would it correspond to "Face set" values if I were to use an exodusii file? Thanks, Justin On Wed, Aug 26, 2015 at 9:17 AM, David Ham <David.Ham@imperial.ac.uk> wrote:
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@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@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@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@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@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@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@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@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@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@imperial.ac.uk >> >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> > >> > >> > _______________________________________________ >> > firedrake mailing list >> > firedrake@imperial.ac.uk >> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> > >> >> _______________________________________________ >> firedrake mailing list >> firedrake@imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > > _______________________________________________ > firedrake mailing list > firedrake@imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 28 Aug 2015, at 00:38, Justin Chang <jychang48@gmail.com> wrote:
So if I use a GMSH file, the third argument in DirichletBC(...) would correspond to whatever the Physical IDS are? Would it correspond to "Face set" values if I were to use an exodusii file?
Yes, that's right. This is meant to work, but I think we have tested the exodus reader less extensively than the gmsh one, so please do complain if things don't work. Cheers, Lawrence
participants (3)
-
David Ham
-
Justin Chang
-
Lawrence Mitchell