Spatially dependent mannings coefficient fails with Nan Flux

[jd-moulton]

We create a mesh with two face sets on the surface to define regions that are then used to define spatially dependent mannings coefficient. This is based on the Borden test case, with a simple rainfall event. The simulation runs if mannings coefficient is defined in the usual way as homogenous over the whole surface. We also tested the regions defined by the two face sets to set a spatially dependent rainfall (different constant over each region). So we are reasonably confident that the face sets themselves, and resulting regions, are well defined.

An input file and mesh file are on Google Drive

Note: I understand this is a crappy location, but issues don't allow these files to be attached (hard to believe I know).

[hydro-geomorph-zhang]

We tested the Borden case simulation with three subsurface layers with different permeability. The first run showed Nan fluxes with a critical depth boundary condition. Then, David suggested making sure a non-zero value for the epsilon coefficient to avoid a zero slope. It did solve the nan flux problem.

The second run was to add different manning's coefficients on the land surface mesh to represent surface heterogeneity. I made up a simple land cover data with only two land cover types to generate two different regions (or facesets). I tried to apply two different rainfall on the two regions with homogeneous manning's N. There is no problem in the simulation, which means that the definition of the two surface regions is correct.

However, after I assigned two different manning's N values on the two regions, the simulation showed an error message "Overland Conductivity Evaluator: Manning coeficient has at least one value that is non-positive. Perhaps you forgot to set the "boundary_face" component?". Apparently, we added boundary_face when we define each manning's N value. Do you think there is anything we missed?

[ecoon]

To attach xml, just change the name to "filename.xml.txt". You can probably do the same with exo files, but that's a little trickier maybe.

[ecoon]

So the problem is cells are 0, not boundary faces. Are you sure that your regions cover the surface?

[jd-moulton]

It is a good question, our other test could have some missing cells from the source and that probably wouldn't matter. I can't do this until tomorrow or Friday, but if you're in the middle of testing could you define the second region as the subtraction of the whole surface with the first set. That would give a complete covering but may still show that there's a problem?

[ecoon]

Yes. We may want to consider having CompositeVectorSpace use MeshPartition instead of MeshFunction. This would confirm that all entities are touched. This was specifically not done on purpose, because it requires uniqueness -- you can't just do something like:

• mesh_function-region1 = computation domain
• mesh_function-region2 = my subset of the domain

(e.g. overwrite the first region with the second region).

But maybe this is less of an issue now -- we primarily used that before RegionLogical was added so that subtraction was doable. So probably we should make an Amanzi ticket to refactor CompositeVectorFunction to use MeshPartition.

[hydro-geomorph-zhang]

I just checked the mesh. Some triangles have one vertex within the domain, but two vertices are out of the domain. So, some centroids are out of the boundary already.

[ecoon]

@jd-moulton see amanzi/amanzi#610

[ecoon]

Yes, this was the problem.

If I add the following regions:

    <ParameterList name="LC_1_PLUS_2" type="ParameterList">
      <ParameterList name="region: logical" type="ParameterList">
        <Parameter name="operation" type="string" value="union" />
        <Parameter name="regions" type="Array(string)" value="{LC_1,LC_2}" />
      </ParameterList>
    </ParameterList>

    <ParameterList name="LC_3" type="ParameterList">
      <ParameterList name="region: logical" type="ParameterList">
        <Parameter name="operation" type="string" value="subtract" />
        <Parameter name="regions" type="Array(string)" value="{computational domain, LC_1_PLUS_2}" />
      </ParameterList>
    </ParameterList>

and then the following Manning coefficient definition:

      <ParameterList name="surface-manning_coefficient" type="ParameterList">

        <Parameter name="field evaluator type" type="string" value="independent variable" />
        <Parameter name="constant in time" type="bool" value="true" />
        <ParameterList name="function" type="ParameterList">

          <ParameterList name="LC_type_1" type="ParameterList">
            <Parameter name="region" type="string" value="LC_1" />
            <Parameter name="components" type="Array(string)" value="{cell,boundary_face}" />
            <ParameterList name="function" type="ParameterList">
              <ParameterList name="function-constant" type="ParameterList">
                <Parameter name="value" type="double" value="0.049" />
              </ParameterList>
            </ParameterList>
          </ParameterList>

          <ParameterList name="LC_type_2" type="ParameterList">
            <Parameter name="region" type="string" value="LC_2" />
            <Parameter name="components" type="Array(string)" value="{cell,boundary_face}" />
            <ParameterList name="function" type="ParameterList">
              <ParameterList name="function-constant" type="ParameterList">
                <Parameter name="value" type="double" value="0.017568" />
              </ParameterList>
            </ParameterList>
          </ParameterList>

          <ParameterList name="LC_type_3" type="ParameterList">
            <Parameter name="region" type="string" value="LC_3" />
            <Parameter name="components" type="Array(string)" value="{cell,boundary_face}" />
            <ParameterList name="function" type="ParameterList">
              <ParameterList name="function-constant" type="ParameterList">
                <Parameter name="value" type="double" value="0.1" />
              </ParameterList>
            </ParameterList>
          </ParameterList>
        </ParameterList>

      </ParameterList>

Then the run starts to run and gets through a few cycles. The image of the Manning coefficient is:

It looks to me like your LC_2 region is empty when resolved on the surface mesh -- did you get the right label? Or is that region defined correctly in the mesh? The exodus file seems to think the region with label 30 has 1575 faces attached to it, but none of them show up in this plot (note the min Manning coef of 0.049, not 0.017568 that I expected).

[jd-moulton]

We'll look more into the mesh. Since Pin was reporting the same problem, it may be worth asking him for his reproducer. It would be good to know if he has the same problem.

[jd-moulton]

Here's a corrected input file, that shows the problem is slightly different. The error I get is

State              |    initializing fields through initial conditions...
terminate called after throwing an instance of 'Errors::Message'
  what():  Overland Conductivity Evaluator: Manning coeficient has at least one value that is non-positive.  Perhaps you forgot to set the "boundary_face" component?

The image below shows that even though the presence of LC_3 allows the simulation to run, no cells are actually set to the value of LC_3. So not sure how / why it is helping.

borden-mix-ethan.xml.txt

[hydro-geomorph-zhang]

Not sure if that would be helpful, but we did find a cell that is out of the domain (see the screenshot attached here). However, it looks like the elevation value, etc was assigned to this cell.

[lipnikov]

FYI, Amanzi has a few tools that helped me tremendously in code debugging. One is the i/o capability State::WriteStateStatistics(S, vo). The other one is mesh region analysis, see Amanzi::InputAnalysis class. You may consider using them in ATS.

[jd-moulton]

With respect to Yu's observation about the data that was used by TINerator to generate the mesh, it appears that TINerator did assign those cells to the appropriate face sets. So I don't think that this is an issue here.

[ecoon]

Ok, good, David's new example makes this clearer. When I run this, with LC_3 included, all cells are either LC_1 or LC_2 value. Without LC_3 included, however, there is at least one boundary face that is 0. So the logic there is off somehow.

[ecoon]

Ok, I'm about halfway through fixing this and realizing that I don't quite trust doing this right now as a minor patch.

The problem is that we currently expect to be able to calculate both Manning's n and slope on both cells and boundary faces.

First, this is an ATS problem, not an Amanzi problem. By this I mean, Amanzi's CompositeVectorFunction::Compute() is correctly implemented. If you try to compute a boundary face value, it uses the corresponding face. For geometric regions, this is the "right" thing to do. For labeled set regions, if those regions are labeled faces (e.g. "surface outlet" region) then this is the right thing to do. But what we want to do in this case is to use not the face with the labeled set but the internal cell. So this is a "special case" and not the right thing to do generically.

Therefore, the better solution is to not use "slope" or "Manning coefficient" vectors from the boundary face, but instead calculate those quantities only on cells and take the boundary face's internal cell's value.

There are some significant structural changes that have to happen to make that work, and I think I have those partially done now. Effectively the overland conductivity evaluator(s) have to become smarter to take the cell value instead of just looping over all components and then all entities in that component and applying the model.

The problem is that this change may affect a bunch of different problems:

• snow redistribution model
• Manning model for freezing case
• Manning subgrid model

I don't have the time to test this sufficiently by tomorrow, as I'm not entirely sure we have the testing in place for any of these. I'd feel better doing this on a branch and including it in next release.

Can you use the following workaround for your demo? Since the boundary faces are all the LC_1 value (with the exception of the outlet) can you use something like:

      <ParameterList name="surface-manning_coefficient" type="ParameterList">

        <Parameter name="field evaluator type" type="string" value="independent variable" />
        <Parameter name="constant in time" type="bool" value="true" />
        <ParameterList name="function" type="ParameterList">

          <ParameterList name="LC_type_1" type="ParameterList">
            <Parameter name="region" type="string" value="LC_1" />
            <Parameter name="components" type="Array(string)" value="{cell}" />
            <ParameterList name="function" type="ParameterList">
              <ParameterList name="function-constant" type="ParameterList">
                <Parameter name="value" type="double" value="0.049" />
              </ParameterList>
            </ParameterList>
          </ParameterList>

          <ParameterList name="LC_type_2" type="ParameterList">
            <Parameter name="region" type="string" value="LC_2" />
            <Parameter name="components" type="Array(string)" value="{cell}" />
            <ParameterList name="function" type="ParameterList">
              <ParameterList name="function-constant" type="ParameterList">
                <Parameter name="value" type="double" value="0.017568" />
              </ParameterList>
            </ParameterList>
          </ParameterList>

          <ParameterList name="LC_type_1_BFs" type="ParameterList">
            <Parameter name="region" type="string" value="surface domain" />
            <Parameter name="components" type="Array(string)" value="{boundary_face}" />
            <ParameterList name="function" type="ParameterList">
              <ParameterList name="function-constant" type="ParameterList">
                <Parameter name="value" type="double" value="0.017568" />
              </ParameterList>
            </ParameterList>
          </ParameterList>

          <ParameterList name="LC_type_2_BFs" type="ParameterList">
            <Parameter name="region" type="string" value="stream outlet" />
            <Parameter name="components" type="Array(string)" value="{boundary_face}" />
            <ParameterList name="function" type="ParameterList">
              <ParameterList name="function-constant" type="ParameterList">
                <Parameter name="value" type="double" value="0.049" />
              </ParameterList>
            </ParameterList>
          </ParameterList>
        </ParameterList>
      </ParameterList>

edit: I've checked and this works for me as a workaround.

[jd-moulton]

For the short course I think we can skip discussing it in the regular session, and have the workaround in an advanced example that participants can play with if they like. So this takes the pressure of having a patch for the short course.

But it would be nice to have this patched on 1.2 eventually, since we would like people to use this branch for upcoming science campaigns. My hope would be that we use 1.2.1 for the short course and provide 1.2.2 to address any issues that crop up during or shortly after the short course - and promote use of 1.2.2 for science campaigns. Would that work?

[ecoon]

Sure, we can try that.

To be clear, I recognize that we need to fix this.

That said -- I have yet to see a run where Manning's n matters. Over daily time scales, Mannings coefficient of < 1 is all just "faster than subsurface." I don't think Manning's n variation affects the results at all until someone does a run that has at least one of the following changes:

1. looks at discharge over time scales << 1 day
2. looks at spatial scales >> 10km.

Ahmad and I have done lots of checking on this on lots of different runs, and Pin has done some checking on this for watershed runs -- it just doesn't change the answer at all for the types of runs we usually do.

[ecoon]

FYI, Amanzi has a few tools that helped me tremendously in code debugging. One is the i/o capability State::WriteStateStatistics(S, vo). The other one is mesh region analysis, see Amanzi::InputAnalysis class. You may consider using them in ATS.

We do use WriteStateStatistics() now, but I'm not familiar with InputAnalysis. Can you give a short description of this, or is it documented somewhere? I don't see any description of what it is doing in the code.

[ecoon]

working on this now on branch ecoon/ism_refactor

[ecoon]

Can @jd-moulton re-check this on Amanzi branch master and ATS branch ecoon/ism_refactor? I hope it should work now.

[ecoon]

that branch is merged into master now, so @jd-moulton and @pinshuai this could be checked on master now. I believe it should work -- feel free to reopen this ticket if you have more trouble.