Prescribe pressure at a corner for Darcy's equation
Hi all, I am trying to simulate the 2D benchmark Elder problem, which is density-driven fluid flow in porous media. The problem is essentially Darcy's equation where no flow boundary conditions are enforced on all four sides (let's suppose the computational domain is 600m by 150m). To ensure uniqueness, pressure is prescribed at the top left and top right corners. Suppose I want to use the RT0 formulation. If my pressure is non-zero, how would I prescribe if I were to use one of the built-in Mesh functions (as opposed to GMSH)? I know for RT0, the pressure bc is prescribed weakly but how do I literally only have the two element faces that connect each of the top left and top right corners apply said boundary condition? Thanks, Justin
Hi Justin, On 01/02/17 09:53, Justin Chang wrote:
Hi all,
I am trying to simulate the 2D benchmark Elder problem, which is density-driven fluid flow in porous media. The problem is essentially Darcy's equation where no flow boundary conditions are enforced on all four sides (let's suppose the computational domain is 600m by 150m). To ensure uniqueness, pressure is prescribed at the top left and top right corners.
Suppose I want to use the RT0 formulation. If my pressure is non-zero, how would I prescribe if I were to use one of the built-in Mesh functions (as opposed to GMSH)? I know for RT0, the pressure bc is prescribed weakly but how do I literally only have the two element faces that connect each of the top left and top right corners apply said boundary condition?
What you need to do is, before doing anything with your mesh, add additional boundary marker labels that correspond to the appropriate facets. You do this half through plex: from firedrake import * mesh = UnitSquareMesh(...) plex = mesh._plex coords = plex.getCoordinates() coord_sec = plex.getCoordinateSection() if plex.getStratumSize("boundary_faces", 1) > 0: faces = plex.getStratumIS("boundary_faces", 1).getIndices() for face in faces: face_coords = plex.vecGetClosure(coord_sec, coords, face) # now face_coords is a numpy array of the vertices of # the face (XY XY) if face_coords are in corner face: plex.setLabelValue("boundary_ids", face, ID) V = FunctionSpace(mesh, "RT", 1) ... now you can use ID just as you do the built in mesh ids. (untested, but this is general idea). If you look in utility_meshes.py at (say) RectangleMesh, you can see how we mark the full faces of the mesh and copy that logic. Lawrence
Thanks Lawrence! I will give this (as well as the other mail about projecting data) a try shortly. On Wed, Feb 1, 2017 at 4:02 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
Hi Justin,
On 01/02/17 09:53, Justin Chang wrote:
Hi all,
I am trying to simulate the 2D benchmark Elder problem, which is density-driven fluid flow in porous media. The problem is essentially Darcy's equation where no flow boundary conditions are enforced on all four sides (let's suppose the computational domain is 600m by 150m). To ensure uniqueness, pressure is prescribed at the top left and top right corners.
Suppose I want to use the RT0 formulation. If my pressure is non-zero, how would I prescribe if I were to use one of the built-in Mesh functions (as opposed to GMSH)? I know for RT0, the pressure bc is prescribed weakly but how do I literally only have the two element faces that connect each of the top left and top right corners apply said boundary condition?
What you need to do is, before doing anything with your mesh, add additional boundary marker labels that correspond to the appropriate facets.
You do this half through plex:
from firedrake import *
mesh = UnitSquareMesh(...)
plex = mesh._plex
coords = plex.getCoordinates() coord_sec = plex.getCoordinateSection()
if plex.getStratumSize("boundary_faces", 1) > 0: faces = plex.getStratumIS("boundary_faces", 1).getIndices() for face in faces: face_coords = plex.vecGetClosure(coord_sec, coords, face) # now face_coords is a numpy array of the vertices of # the face (XY XY) if face_coords are in corner face: plex.setLabelValue("boundary_ids", face, ID)
V = FunctionSpace(mesh, "RT", 1) ...
now you can use ID just as you do the built in mesh ids.
(untested, but this is general idea).
If you look in utility_meshes.py at (say) RectangleMesh, you can see how we mark the full faces of the mesh and copy that logic.
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Lawrence, I tried following the syntax of what I see in utility_meshes.py but for some reason it seems the boundary_ids are not registering. Attached is my code. Run it as: python 2D_Darcy_Homo_Iso.py 600 150 0 Thanks, Justin PS - In the attached code, I could also pass in either -flow_ksp_monitor or -tran_ksp_monitor to view the ksp monitors for the Darcy or transport equations, respectively, but if I try to pass in both, i get this error: Traceback (most recent call last): File "2D_Darcy_Homo_Iso.py", line 197, in <module> solver_parameters=flow_parameters) File "/Users/jychang48/Software/firedrake/src/firedrake/firedrake/variational_solver.py", line 262, in __init__ super(LinearVariationalSolver, self).__init__(*args, **kwargs) File "/Users/jychang48/Software/firedrake/src/firedrake/firedrake/variational_solver.py", line 165, in __init__ self.set_from_options(self.snes) File "/Users/jychang48/Software/firedrake/src/firedrake/firedrake/solving_utils.py", line 102, in set_from_options petsc_obj.setFromOptions() File "PETSc/SNES.pyx", line 112, in petsc4py.PETSc.SNES.setFromOptions (src/petsc4py.PETSc.c:162247) petsc4py.PETSc.Error: error code 56 [0] SNESSetFromOptions() line 983 in /private/var/folders/92/1kh0g4kn2z50fnwsmf8s269w0000gn/T/pip-S9HxIn-build/src/snes/interface/snes.c [0] KSPSetFromOptions() line 499 in /private/var/folders/92/1kh0g4kn2z50fnwsmf8s269w0000gn/T/pip-S9HxIn-build/src/ksp/ksp/interface/itcl.c [0] KSPMonitorSetFromOptions() line 321 in /private/var/folders/92/1kh0g4kn2z50fnwsmf8s269w0000gn/T/pip-S9HxIn-build/src/ksp/ksp/interface/itcl.c [0] PetscOptionsGetViewer() line 241 in /private/var/folders/92/1kh0g4kn2z50fnwsmf8s269w0000gn/T/pip-S9HxIn-build/src/sys/classes/viewer/interface/viewreg.c [0] No support for this operation for this object type [0] Unsupported viewer -tran_ksp_monitor I can't seem to pinpoint what might be causing this. On Wed, Feb 1, 2017 at 4:09 AM, Justin Chang <jychang48@gmail.com> wrote:
Thanks Lawrence! I will give this (as well as the other mail about projecting data) a try shortly.
On Wed, Feb 1, 2017 at 4:02 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
Hi Justin,
On 01/02/17 09:53, Justin Chang wrote:
Hi all,
I am trying to simulate the 2D benchmark Elder problem, which is density-driven fluid flow in porous media. The problem is essentially Darcy's equation where no flow boundary conditions are enforced on all four sides (let's suppose the computational domain is 600m by 150m). To ensure uniqueness, pressure is prescribed at the top left and top right corners.
Suppose I want to use the RT0 formulation. If my pressure is non-zero, how would I prescribe if I were to use one of the built-in Mesh functions (as opposed to GMSH)? I know for RT0, the pressure bc is prescribed weakly but how do I literally only have the two element faces that connect each of the top left and top right corners apply said boundary condition?
What you need to do is, before doing anything with your mesh, add additional boundary marker labels that correspond to the appropriate facets.
You do this half through plex:
from firedrake import *
mesh = UnitSquareMesh(...)
plex = mesh._plex
coords = plex.getCoordinates() coord_sec = plex.getCoordinateSection()
if plex.getStratumSize("boundary_faces", 1) > 0: faces = plex.getStratumIS("boundary_faces", 1).getIndices() for face in faces: face_coords = plex.vecGetClosure(coord_sec, coords, face) # now face_coords is a numpy array of the vertices of # the face (XY XY) if face_coords are in corner face: plex.setLabelValue("boundary_ids", face, ID)
V = FunctionSpace(mesh, "RT", 1) ...
now you can use ID just as you do the built in mesh ids.
(untested, but this is general idea).
If you look in utility_meshes.py at (say) RectangleMesh, you can see how we mark the full faces of the mesh and copy that logic.
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Sorry, forgot to mention this: For ID 5, I want this applied to all faces that have y == Ly and (150 <= x <= 450) For ID 6, I want this applied to only the top left and top right most element faces (i.e., the horizontal and vertical components). On Mon, Feb 6, 2017 at 10:11 PM, Justin Chang <jychang48@gmail.com> wrote:
Lawrence,
I tried following the syntax of what I see in utility_meshes.py but for some reason it seems the boundary_ids are not registering.
Attached is my code. Run it as:
python 2D_Darcy_Homo_Iso.py 600 150 0
Thanks, Justin
PS - In the attached code, I could also pass in either -flow_ksp_monitor or -tran_ksp_monitor to view the ksp monitors for the Darcy or transport equations, respectively, but if I try to pass in both, i get this error:
Traceback (most recent call last): File "2D_Darcy_Homo_Iso.py", line 197, in <module> solver_parameters=flow_parameters) File "/Users/jychang48/Software/firedrake/src/firedrake/ firedrake/variational_solver.py", line 262, in __init__ super(LinearVariationalSolver, self).__init__(*args, **kwargs) File "/Users/jychang48/Software/firedrake/src/firedrake/ firedrake/variational_solver.py", line 165, in __init__ self.set_from_options(self.snes) File "/Users/jychang48/Software/firedrake/src/firedrake/firedrake/solving_utils.py", line 102, in set_from_options petsc_obj.setFromOptions() File "PETSc/SNES.pyx", line 112, in petsc4py.PETSc.SNES.setFromOptions (src/petsc4py.PETSc.c:162247) petsc4py.PETSc.Error: error code 56 [0] SNESSetFromOptions() line 983 in /private/var/folders/92/ 1kh0g4kn2z50fnwsmf8s269w0000gn/T/pip-S9HxIn-build/src/snes/ interface/snes.c [0] KSPSetFromOptions() line 499 in /private/var/folders/92/ 1kh0g4kn2z50fnwsmf8s269w0000gn/T/pip-S9HxIn-build/src/ksp/ ksp/interface/itcl.c [0] KSPMonitorSetFromOptions() line 321 in /private/var/folders/92/ 1kh0g4kn2z50fnwsmf8s269w0000gn/T/pip-S9HxIn-build/src/ksp/ ksp/interface/itcl.c [0] PetscOptionsGetViewer() line 241 in /private/var/folders/92/ 1kh0g4kn2z50fnwsmf8s269w0000gn/T/pip-S9HxIn-build/src/sys/ classes/viewer/interface/viewreg.c [0] No support for this operation for this object type [0] Unsupported viewer -tran_ksp_monitor
I can't seem to pinpoint what might be causing this.
On Wed, Feb 1, 2017 at 4:09 AM, Justin Chang <jychang48@gmail.com> wrote:
Thanks Lawrence! I will give this (as well as the other mail about projecting data) a try shortly.
On Wed, Feb 1, 2017 at 4:02 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
Hi Justin,
On 01/02/17 09:53, Justin Chang wrote:
Hi all,
I am trying to simulate the 2D benchmark Elder problem, which is density-driven fluid flow in porous media. The problem is essentially Darcy's equation where no flow boundary conditions are enforced on all four sides (let's suppose the computational domain is 600m by 150m). To ensure uniqueness, pressure is prescribed at the top left and top right corners.
Suppose I want to use the RT0 formulation. If my pressure is non-zero, how would I prescribe if I were to use one of the built-in Mesh functions (as opposed to GMSH)? I know for RT0, the pressure bc is prescribed weakly but how do I literally only have the two element faces that connect each of the top left and top right corners apply said boundary condition?
What you need to do is, before doing anything with your mesh, add additional boundary marker labels that correspond to the appropriate facets.
You do this half through plex:
from firedrake import *
mesh = UnitSquareMesh(...)
plex = mesh._plex
coords = plex.getCoordinates() coord_sec = plex.getCoordinateSection()
if plex.getStratumSize("boundary_faces", 1) > 0: faces = plex.getStratumIS("boundary_faces", 1).getIndices() for face in faces: face_coords = plex.vecGetClosure(coord_sec, coords, face) # now face_coords is a numpy array of the vertices of # the face (XY XY) if face_coords are in corner face: plex.setLabelValue("boundary_ids", face, ID)
V = FunctionSpace(mesh, "RT", 1) ...
now you can use ID just as you do the built in mesh ids.
(untested, but this is general idea).
If you look in utility_meshes.py at (say) RectangleMesh, you can see how we mark the full faces of the mesh and copy that logic.
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hi Justin,
On 7 Feb 2017, at 04:13, Justin Chang <jychang48@gmail.com> wrote:
Sorry, forgot to mention this:
For ID 5, I want this applied to all faces that have y == Ly and (150 <= x <= 450)
For ID 6, I want this applied to only the top left and top right most element faces (i.e., the horizontal and vertical components).
On Mon, Feb 6, 2017 at 10:11 PM, Justin Chang <jychang48@gmail.com> wrote: Lawrence,
I tried following the syntax of what I see in utility_meshes.py but for some reason it seems the boundary_ids are not registering
I had a look, I think the problem is like so. Our code assumes that each mesh facet only ever has one boundary id. So what's happening is that the plex is correctly getting the new boundary ids, but when we look them up the ones that were already created appear "first" and override your new ones. What you can do is clear all the boundary id labels and then explicitly put in all the IDs you want (in such a way that they don't overlap). I think the code could probably support multiple labels on each facet, but would need to unpick it a bit: can you raise an issue?
Attached is my code. Run it as:
python 2D_Darcy_Homo_Iso.py 600 150 0
Thanks, Justin
PS - In the attached code, I could also pass in either -flow_ksp_monitor or -tran_ksp_monitor to view the ksp monitors for the Darcy or transport equations, respectively, but if I try to pass in both, i get this error:
This one is easy (!). We inspect the PETSc options database like so: from petsc4py import PETSc PETSc.Options(prefix).getAll() If I have options -foo_bar -baz_bar then if I do: PETSc.Options("foo_").getAll() => {"bar": "-baz_bar"} I will push a fix for this one. Lawrence
On 7 Feb 2017, at 09:03, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:
PETSc.Options("foo_").getAll()
=> {"bar": "-baz_bar"}
I will push a fix for this one.
Please try the branch proposed here: https://github.com/firedrakeproject/firedrake/pull/1004 Lawrence
Hi again, So this may or may not be related to the original issue. I was able to define corner IDs needed for my pressure BCs but I cannot seem to enforce both no-flow and atmospheric pressure simultaneously. Enforcing no-flow seems to negate the weakly prescribed pressure flux and my pressure solution is complete wack. On the otherhand, not enforcing no flow gives me a pressure profile I expect, but I no longer have a no flow domain. Attached is my code, the custom boundary IDs begin at line 39, the weak formulations begin at line 164. Run the code as: python 2D_Darcy_Homo_Iso.py 600 150 2 The code right now will generate a nonsensical pressure solution because I have no-flow enforced on all boundaries (including the corners). Any ideas on what to do here? Thanks, Justin On Tue, Feb 7, 2017 at 4:55 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
On 7 Feb 2017, at 09:03, Lawrence Mitchell <lawrence.mitchell@imperial. ac.uk> wrote:
PETSc.Options("foo_").getAll()
=> {"bar": "-baz_bar"}
I will push a fix for this one.
Please try the branch proposed here:
https://github.com/firedrakeproject/firedrake/pull/1004
Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (2)
-
Justin Chang
-
Lawrence Mitchell