Observation region using WW

[Chuyang Liu]

Hi there,

Can you guide me through the process of designating specific elements in the mesh as regions or sidesets using the Watershed Workflow? For instance, I have a monitoring well situated within the simulation domain. I want to focus on the regions that encompass the well's location – spanning from the ground surface to the bottom of the simulation domain or a designated mesh element at a specific depth. I'm uncertain about how to translate the geographic coordinates of the well into corresponding sidesets within the mesh file in a straightforward way.

Thanks,

Chuyang

[Coon, Ethan]

Good question.  We don’t have a good example of doing this, but it is fairly straightforward conceptually, and a common enough thing to need!

 

Also, you can do this in two ways – either in WW to set up a region, or in ATS directly.

 

In ATS directly:

So ATS regions can be defined using geometry e.g. boxes and lines and more.  I think probably the easiest bet is to define your well as a “region: line segment”:

 

https://amanzi.github.io/ats/input_spec/ATSNativeSpec_dev.html#line-segment

 

which was intended for wells.  Just set the top and bottom points and you’ll get all cells that are touching that line.

 

In WW as a LabeledSet:

 

A region is created using ww.mesh.LabeledSet.  The mesh has a list of them, so you just want to create the LabeledSet and append it to that list, e.g, assuming you have a mesh M, you might:

 

new_id = M.next_available_labeled_setid()

M.labeled_sets.append(LabeledSet(‘my_well’, new_id, ‘CELL’, [list_of_IDs])

 

So there are two ways you can get the list of IDs.

 

First, let’s assume you’ve identified the surface cell that includes the horizontal coordinates.  You should do this through something like (noting that I am making this up as I go, so there are probably bugs here):

 

my_well_xy = [x,y]

my_well_cell = -1

for ic, cell in enumerate(m2.conn):

        points = np.array([m2.coords[i] for i in cell])

        if shapely.geometry.Polygon(points).contains(my_well_xy):

            my_well_cell = ic

            break

 

So now we have the 2D cell. 

 

If you want the entire column of IDs, then the easiest thing is to put a LabeledSet on the 2D mesh, using the option “to_extrude=True” when constructing the LabeledSet.  This tells the code that, when you extrude from 2D to 3D, to apply this labeled set on all cells that are extruded.  (If this is false, it only keeps the top face of the column).

 

If, instead, you want a subset of cells in the column, then you need to add the LabeledSet to the 3D mesh, and you need the cell IDs in the vertical column.  Cells in a column are ordered continuously, so the formula is:

 

3D_cell_id = z_cell + 2D_cell_id * ncells_extruded

 

Where z_cell is the 0 (top) to ncells_extruded -1 (bottom) cell in the vertical.

 

Given a depth z, you would have to go through your dz used in extrusion to figure out which cell it was (z_cell).  Once you have that list, you can make a LabeledSet.

 

Ethan

 

 

 

 

 

 

 

--
You received this message because you are subscribed to the Google Groups "Amanzi-ATS Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ats-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/ats-users/5be8ff18-d7ea-479e-8fb9-111f008eef5cn%40googlegroups.com.

[Chuyang Liu]

Thank you very much, Ethan. With your suggestions, I've successfully created and tested a Python function. By the way, I have another question regarding the mesh. Is it feasible to have an uneven ground surface with a flat surface at the bottom?

[Ethan Coon]

Yes.  When you use the extrude function, most workflows pass in the option “CONSTANT” as the layer_type (as in constant dz).  If you look at the documentation for the extrusion function, you’ll see that there are a few choices there, including one to set a flat bottom layer at a constant elevation, which will set dz as needed to hit that elevation.  Alternatively, you can provide the dz for every stack of nodes and take as much control as you want.

Ethan

To view this discussion on the web visit https://groups.google.com/d/msgid/ats-users/3f972023-6ad4-48bd-89f9-d8ad38a85ea1n%40googlegroups.com.

--

-------------------------------------------------------------------
Ethan Coon
917-969-6831
https://www.ornl.gov/staff-profile/ethan-t-coon
-------------------------------------------------------------------