Problems in Erobility strcture with 3D displacements

[Picus-canus]

I am currently working on a 3D model of uwgeodynamics, whose surface processes need to be coupled to badlands, and I want to give the region a different K value.

When running the code, I noticed that I would be placements at function apply_XY_displacements

This is because it is trying to assign a two-dimensional array to a one-dimensional array.

Combined with the code before and after this function, I think this paragraph is calculating properties for the new node in TIN, so I changed Tevals to a 3D array(line 911):
Tevals= numpy.zeros((len(zvals), ids.shape[1], lay))

Change the Tevals assignment to(line 918):
Tevals[:, k, :] = mTe[ids[:, k], :]

Average instead(line 934):
Teavg[:, k] = numpy.average(Tevals[:, :, k], weights=weights, axis=1)

It is now ready to run successfully from apply_XY_displacements
However, in flow.erodibility[1GIDs] in checkpoints, Th and Ke values contain addPts, resulting in the following errors:

I think this may have to do with the assignment to flow.erodibility in the _rebuild_mesh because mapero.erodibility is still a one-dimensional array without addPts and the assignment is not the same as the module above.

But I don't know much about the logic of this block and would like some help.
Thank you for your time. I appreciate any help and look forward to the reply.

[tristan-salles]

Hi @Picus-canus, from memory the coupled badlands/UDW does not allow for variable erodibility map... which might be the reason why you are having the issue.

After your modifications, I think Kevals or Keavg that you are updating based on the new point should be assigned to the variable flow.erodibility which is initially set based on self.mapero.erodibility. I haven't tested it but this might help to debug the code...

[Picus-canus]

Thank you four your reply, @tristan-salles !

Yes, I noticed that the XML input description mentioned that 3D displacements and variable erodibility map would have issues - after I had actually tried to make them.

Btw, I think this may be the last step of the edit of code, for it is the time to checkpoint and about to export h5py files. I had tested the array size of self.mapero.erodibility, and it did not update the addPTs of the new ids index, taking the example in my issue graph2, self.mapero.erodibility has the shape (31966,), This also leads to the shape of self.flow.erodibility (31966,).

I think using the eroMesh.eroMesh function will work because mapero has used this function before (_build_mesh-buildMesh.construct_mesh-_build_strateroMesh-eroMesh.eroMesh), and it updates xyTIN.

Hope this function can be debugged successfully! Thank you again and looking forward to your reply.

[tristan-salles]

Hi again,

• Yes this is indeed happening at the last step prior to h5py files export.
• So you have the stratigraphy turned on as well? I am asking because eroMesh.eroMesh defines the underlying stratigraphy.

Maybe what also needs to be changed is at line 935 in the apply_XY_displacements function:

Keavg[:, k] = numpy.average(Keavg[:, :, k], weights=weights, axis=1)

return (
            newTIN,
            newelev,
            newcum,
            newhcum,
            newfcum,
            newwcum,
            newcumf,
            newscum,
            newKe,
            newTe,
        )

All the return variables here should have the same dimensions.

Then in model.py, the newKe is set to self.mapero.Ke and you might need to force flow.erodibility to be set to the value of the mapero.erodibility following eroMesh.py (line 80 to 90):

            # Get erodibility from erosive layer thicknesses
            self.erodibility = numpy.zeros(nbPts)
            for k in range(self.layNb):
                existIDs = numpy.where(
                    numpy.logical_and(
                        self.thickness[:, k] > 0.0, self.erodibility[:] == 0.0
                    )
                )[0]
                self.erodibility[existIDs] = self.Ke[existIDs, k]
                if len(numpy.where(self.erodibility == 0)[0]) == 0:
                    break 

Basically self.flow.erodibility is the top layer or non-null deposit values of self.mapero.erodibility.

[Picus-canus]

Thanks agian, @tristan-salles ! I give myself a weekend off so I'm only replying now :)

Yes, I modified Ke in the same way, because Th and Ke are both 2D arrays on apply_XY_displacements, so Kevals was also modified as 3D arrays and Keavg was processed in the same way as Teavg. Only in this way can it successfully run apply_XY_displacements.

# 912 
Kevals = numpy.zeros((len(zvals), ids.shape[1], lay))
# 919 
Kevals[:, k, :] = mTe[ids[:, k], :]
# 935 
Keavg[:, k] = numpy.average(Kevals[:, :, k], weights=weights, axis=1)

I think the code you pointed out to me is really the one I need right now!
But eroMesh is like a function that is designed to be used when building a mesh, and its input is a bit complicated for me, so I decide to write a new function.

# in eroMesh.py
def update_erodibility(self, erodibility, Ke, thickness):
    self.erodibility = numpy.zeros(len(Ke))
    for k in range(self.Ke.shape[1]):
        existIDs = numpy.where(
            numpy.logical_and(
                self.thickness[:, k] > 0.0, self.erodibility[:] == 0.0
            )
        )[0]
        self.erodibility[existIDs] = self.Ke[existIDs, k]
        if len(numpy.where(self.erodibility == 0)[0]) == 0:
            break

and also call this function in Model.py line 267:

#line 266
        else:
#line 267-270
            self.mapero.update_erodibility(
                self.mapero.erodibility,
                self.mapero.Ke,
                self.mapero.thickness
            )
#line 271
            self.flow.erodibility = self.mapero.erodibility

The good news is that it can now complete a loop with no errors, but now it has a new problem.
Again in apply_XY_displacements, it seems newKe variable is not created.
When I'm looking at the if condition, this only happens if len(newTIN["vertices"][:, 0]) = len(newXY[:, 0]), but I haven't done debugging yet, and I'm not sure if I can just change it to >=.

As the cycle can be completed, I believe it should be close to success. Thank you again and looking forward to your reply.

[Picus-canus]

I test change
if len(newTIN["vertices"][:, 0]) > len(newXY[:, 0])
to
if len(newTIN["vertices"][:, 0]) >= len(newXY[:, 0])

and it can loop more times, with next problem...

I will track the shape of Ke tomorrow.

Besides, what I'm testing is to add the following to the "examples/rift/rift.xml" file

    <erocoeff>
        <!-- Number of erosion layers. -->
        <erolayers>1</erolayers>
        <!-- The layering is defined from top to bottom, with:
            - either a constant erodibility value for the entire layer or with an erodibility map
            - either a constant thickness for the entire layer or with a thickness map -->
        <!-- Constant erodibility and layer thickness -->
        <erolay>
            <!-- Uniform erodibility value for the considered layer. -->
            <erocst>1.e-6</erocst>
            <!-- Uniform thickness value for the considered layer [m]. -->
            <thcst>1000</thcst>
        </erolay>
    </erocoeff>

I also used my model running the code I programmed, the same problem turned out:

I have no idea why the array of Ke is not the same as the array of elev.