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