[firedrake] Prescribe pressure at a corner for Darcy's equation

[Justin Chang]

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 at 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 at 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 at 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 at imperial.ac.uk
>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>>>
>>>
>>
>
-------------- next part --------------
HTML attachment scrubbed and removed