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