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

[Justin Chang]

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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 2D_Darcy_Homo_Iso.py
Type: text/x-python-script
Size: 6719 bytes
Desc: not available
URL: <http://mailman.ic.ac.uk/pipermail/firedrake/attachments/20170206/44a5dc90/attachment.bin>