Re: [firedrake] Laplace equation on subdomain
Hi Thomasz, You can get the matching CG1 indicator function (to use in the custom BC) using: I_cg = Function(V) par_loop('for (int i=0; i<A.dofs; i++) A[i][0] = fmax(A[i][0], B[0][0]);', dx, {'A' : (I_cg, RW), 'B': (I_f, READ)}) This sets nodes of l_cg to 1 if they touch any cell for which the indicator function is 1. Regards, David On Fri, 15 Apr 2016 at 13:43 Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:
Ok, I extended the BCs. I'm getting indices out of bounds error. I suppose this is because I enforce the BCs in DG0, and then solve in CG space. How can I map it from one to another?
Thanks, Tomasz
Code:
from firedrake import * import numpy as np
nx = 8; ny = 4; Lx = 2.; Ly = 1.
mesh = RectangleMesh( nx, ny, Lx, Ly ) V = FunctionSpace( mesh, "CG", 1 ) V_DG0 = FunctionSpace( mesh, "DG", 0 )
def Heaviside( zero_on_left=True, Lp=0. ): sign = 1. if zero_on_left == False: sign = -1. Hx = Function( V_DG0 ) Hx_expr = Expression( "0.5*(1.0 + copysign( 1.0, sign*(x[0]-Lp) ))", sign=sign, Lp=Lp) Hx.interpolate(Hx_expr) return Hx;
I_f = Heaviside( zero_on_left=False, Lp=1. )
phi = Function( V ) trial = TrialFunction( V ) v = TestFunction( V )
a = dot( grad(trial), grad(v) ) * dx L = Constant(0.) * v * dx
bc = DirichletBC( V, Expression( "sin(pi*x[1]/Ly)", Ly=Ly ), 1 )
class MyBC(DirichletBC): def __init__(self, V, value, markers): # Call superclass init # We provide a dummy subdomain id. super(MyBC, self).__init__(V, value, 0) # Override the "nodes" property which says where the boundary # condition is to be applied. self.nodes = np.unique(np.where(markers.dat.data_ro_with_halos == 0)[0])
# Now I can use MyBC to create a "boundary condition" to zero out all # the nodes that are *not* in the fluid domain: BC_exclude_beyond_fluid = MyBC( V, 0, I_f )
solve( a==L, phi, bcs=[bc, BC_exclude_beyond_fluid] )
Error listing:
Traceback (most recent call last): File "laplace_extended_domain.py", line 51, in <module> solve( a==L, phi, bcs=[bc, BC_exclude_beyond_fluid] ) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/firedrake/solving.py", line 120, in solve _solve_varproblem(*args, **kwargs) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/firedrake/solving.py", line 147, in _solve_varproblem solver.solve() File "<decorator-gen-295>", line 2, in solve File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/profiling.py", line 203, in wrapper return f(*args, **kwargs) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/firedrake/variational_solver.py", line 181, in solve bc.apply(self._problem.u) File "<decorator-gen-297>", line 2, in apply File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/profiling.py", line 203, in wrapper return f(*args, **kwargs) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/firedrake/bcs.py", line 233, in apply r.assign(self.function_arg, subset=self.node_set) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/utils.py", line 64, in __get__ obj.__dict__[self.__name__] = result = self.fget(obj) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/firedrake/bcs.py", line 162, in node_set return op2.Subset(self._function_space.node_set, self.nodes) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/backends.py", line 118, in __call__ return t(*args, **kwargs) File "<decorator-gen-7>", line 2, in __init__ File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/utils.py", line 130, in wrapper return f(*args, **kwargs) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/base.py", line 840, in __init__ (self._indices[0], self._indices[-1], self._superset.total_size)) pyop2.exceptions.SubsetIndexOutOfBounds: Out of bounds indices in Subset construction: [20, 63) not [0, 45)
________________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> Sent: 15 April 2016 11:22 To: firedrake@imperial.ac.uk Subject: Re: [firedrake] Laplace equation on subdomain
On 15 Apr 2016, at 11:08, Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:
Dear Firedrakers,
As an intermediate part of fluid-structure interaction problem, I would like to solve Laplace equation (with zero Neumann boundary conditions) on a subdomain of the full domain. For that reason I define the step function that is 1 in the subdomain of interest and 0 elsewhere. Multiplying LHS of the weak form (with trial function) by the step function causes a solver error. How can I correct it? Code and error listing are below.
The problem is that you've knocked out a part of your operator so that it's:
laplace 0 0 0
Particularly, it's zero on the diagonal outside the subdomain.
To correct this, you'll need to define your own custom boundary condition class that inherits from DirichletBC that identifies the nodes that are outside the subdomain and impose strong conditions there. I think we discussed this when you visited imperial a while ago.
That way you'll end up with:
laplace 0 0 I
with the identity on the diagonal outside the subdomain instead.
As to your comment (how to write 0*v*dx more succinctly) use:
L = Constant(0)*v*dx
Thanks,
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
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? Regards, Tomasz ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of David Ham <David.Ham@imperial.ac.uk> Sent: 15 April 2016 13:17 To: firedrake Subject: Re: [firedrake] Laplace equation on subdomain Hi Thomasz, You can get the matching CG1 indicator function (to use in the custom BC) using: I_cg = Function(V) par_loop('for (int i=0; i<A.dofs; i++) A[i][0] = fmax(A[i][0], B[0][0]);', dx, {'A' : (I_cg, RW), 'B': (I_f, READ)}) This sets nodes of l_cg to 1 if they touch any cell for which the indicator function is 1. Regards, David On Fri, 15 Apr 2016 at 13:43 Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk<mailto:mmtjs@leeds.ac.uk>> wrote: Ok, I extended the BCs. I'm getting indices out of bounds error. I suppose this is because I enforce the BCs in DG0, and then solve in CG space. How can I map it from one to another? Thanks, Tomasz Code: from firedrake import * import numpy as np nx = 8; ny = 4; Lx = 2.; Ly = 1. mesh = RectangleMesh( nx, ny, Lx, Ly ) V = FunctionSpace( mesh, "CG", 1 ) V_DG0 = FunctionSpace( mesh, "DG", 0 ) def Heaviside( zero_on_left=True, Lp=0. ): sign = 1. if zero_on_left == False: sign = -1. Hx = Function( V_DG0 ) Hx_expr = Expression( "0.5*(1.0 + copysign( 1.0, sign*(x[0]-Lp) ))", sign=sign, Lp=Lp) Hx.interpolate(Hx_expr) return Hx; I_f = Heaviside( zero_on_left=False, Lp=1. ) phi = Function( V ) trial = TrialFunction( V ) v = TestFunction( V ) a = dot( grad(trial), grad(v) ) * dx L = Constant(0.) * v * dx bc = DirichletBC( V, Expression( "sin(pi*x[1]/Ly)", Ly=Ly ), 1 ) class MyBC(DirichletBC): def __init__(self, V, value, markers): # Call superclass init # We provide a dummy subdomain id. super(MyBC, self).__init__(V, value, 0) # Override the "nodes" property which says where the boundary # condition is to be applied. self.nodes = np.unique(np.where(markers.dat.data_ro_with_halos == 0)[0]) # Now I can use MyBC to create a "boundary condition" to zero out all # the nodes that are *not* in the fluid domain: BC_exclude_beyond_fluid = MyBC( V, 0, I_f ) solve( a==L, phi, bcs=[bc, BC_exclude_beyond_fluid] ) Error listing: Traceback (most recent call last): File "laplace_extended_domain.py", line 51, in <module> solve( a==L, phi, bcs=[bc, BC_exclude_beyond_fluid] ) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/firedrake/solving.py", line 120, in solve _solve_varproblem(*args, **kwargs) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/firedrake/solving.py", line 147, in _solve_varproblem solver.solve() File "<decorator-gen-295>", line 2, in solve File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/profiling.py", line 203, in wrapper return f(*args, **kwargs) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/firedrake/variational_solver.py", line 181, in solve bc.apply(self._problem.u) File "<decorator-gen-297>", line 2, in apply File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/profiling.py", line 203, in wrapper return f(*args, **kwargs) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/firedrake/bcs.py", line 233, in apply r.assign(self.function_arg, subset=self.node_set) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/utils.py", line 64, in __get__ obj.__dict__[self.__name__] = result = self.fget(obj) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/firedrake/bcs.py", line 162, in node_set return op2.Subset(self._function_space.node_set, self.nodes) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/backends.py", line 118, in __call__ return t(*args, **kwargs) File "<decorator-gen-7>", line 2, in __init__ File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/utils.py", line 130, in wrapper return f(*args, **kwargs) File "/home/tommy/programs/firedrake/local/lib/python2.7/site-packages/pyop2/base.py", line 840, in __init__ (self._indices[0], self._indices[-1], self._superset.total_size)) pyop2.exceptions.SubsetIndexOutOfBounds: Out of bounds indices in Subset construction: [20, 63) not [0, 45) ________________________________________ From: firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk> <firedrake-bounces@imperial.ac.uk<mailto:firedrake-bounces@imperial.ac.uk>> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk<mailto:lawrence.mitchell@imperial.ac.uk>> Sent: 15 April 2016 11:22 To: firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> Subject: Re: [firedrake] Laplace equation on subdomain
On 15 Apr 2016, at 11:08, Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk<mailto:mmtjs@leeds.ac.uk>> wrote:
Dear Firedrakers,
As an intermediate part of fluid-structure interaction problem, I would like to solve Laplace equation (with zero Neumann boundary conditions) on a subdomain of the full domain. For that reason I define the step function that is 1 in the subdomain of interest and 0 elsewhere. Multiplying LHS of the weak form (with trial function) by the step function causes a solver error. How can I correct it? Code and error listing are below.
The problem is that you've knocked out a part of your operator so that it's: laplace 0 0 0 Particularly, it's zero on the diagonal outside the subdomain. To correct this, you'll need to define your own custom boundary condition class that inherits from DirichletBC that identifies the nodes that are outside the subdomain and impose strong conditions there. I think we discussed this when you visited imperial a while ago. That way you'll end up with: laplace 0 0 I with the identity on the diagonal outside the subdomain instead. As to your comment (how to write 0*v*dx more succinctly) use: L = Constant(0)*v*dx Thanks, Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 15 Apr 2016, at 15:54, Tomasz Salwa [RPG] <mmtjs@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
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@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> Sent: 15 April 2016 16:00 To: firedrake@imperial.ac.uk Subject: Re: [firedrake] Laplace equation on subdomain
On 15 Apr 2016, at 15:54, Tomasz Salwa [RPG] <mmtjs@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
participants (3)
-
David Ham
-
Lawrence Mitchell
-
Tomasz Salwa [RPG]