(Click for larger image)
I've got the current C++ topographic steering code converted into Python, which is how I'm doing the testing outside of Child. Before I try the new JP (and other) code, I wanted to make sure I was getting the same results from meander(). At first I'm just capturing the inputs and outputs for each meandering reach in CalcMigration(), but eventually I will add some kind of 'hot link' between the C++ and Python code so that I can continue to run and test transparently in the hybrid system. Once I have the new code working, I'll translate it back into C++ for static linking into Child like the current code.
The first thing I wanted to try was the 'control' case, where the stream flow is linear and constant over time, with no change in (x, y) coordinates. However, when I try running on the smooth dihedral mesh as before, but with the meandering option on, I get the following aberrant behavior. The input fie is here:
The tail end of the terminal output is:
... Writing data for time zero... tOutput::WriteOutput() time=0 Initialization done. 0 0 tOutput::WriteOutput() time=1 1 1 tOutput::WriteOutput() time=2 2 2 tOutput::WriteOutput() time=3 3 3 tOutput::WriteOutput() time=4 4 4 tOutput::WriteOutput() time=5 5 5 tOutput::WriteOutput() time=6 6 6 ../Code/MeshElements/meshElements.cpp:442: failed assertion `varea>0.0' Abort
Here're plots for the first 6 time steps:
(Click for larger image)
I thought that maybe this had something to do with the linearity of the mesh and stream and resulting zero curvature for all meandering nodes, so I tried perturbing the (x, y) coords randomly by (±17, ±20), which eliminated the error but not the removal of interior nodes:
(Click for larger image)
Here's the end result of 1000 steps:
I then tried perturbing the elevation randomly by ±1m:
... Writing data for time zero... tOutput::WriteOutput() time=0 Initialization done. 0 0 tOutput::WriteOutput() time=1 1 1 WARNING-Type 1, from tStreamNet::CalcSlopes....detected a meander node without positive drainage ID= 281, X= 450.333, Y= 180, Z= 1.43522 WARNING-Type 2, from tStreamNet::CalcSlopes....detected a meander node without positive drainage ID= 281, X= 450.333, Y= 180, Z= 1.43522 WARNING-Type 1, from tStreamNet::CalcSlopes....detected a meander node without positive drainage ID= 285, X= 450.333, Y= 740, Z= 6.99306 WARNING-Type 2, from tStreamNet::CalcSlopes....detected a meander node without positive drainage ID= 285, X= 450.333, Y= 740, Z= 6.99306 ... tOutput::WriteOutput() time=13 13 13 WARNING-Type 1, from tStreamNet::CalcSlopes....detected a meander node without positive drainage ID= 89, X= 484.974, Y= 640, Z= 5.5918 WARNING-Type 2, from tStreamNet::CalcSlopes....detected a meander node without positive drainage ID= 89, X= 484.974, Y= 640, Z= 5.5918 WARNING-Type 1, from tStreamNet::CalcSlopes....detected a meander node without positive drainage ID= 89, X= 484.974, Y= 640, Z= 5.5918 WARNING-Type 2, from tStreamNet::CalcSlopes....detected a meander node without positive drainage ID= 89, X= 484.974, Y= 640, Z= 5.5918 ../Code/MeshElements/meshElements.cpp:442: failed assertion `varea>0.0' Abort
(Click for larger image)
(Although I can see why it's desirable for the meandering nodes to have a downward slope, I don't think this is reasonable for a real river, since the bottom depth can change considerably from point to point while the water surface still maintains a downward slope. In a sense, a river is like a lake where every downstream node is an outlet. I may have to do something about this for the MNRR simulation. There's also the issue of width: I don't think it will be sufficient to model the MNRR by a single ideal centerline of nodes. I may need to create a way for the river to span several lateral nodes at each 'point' downstream, especially if we are also going to activate the floodplain option. I'll look into this as it becomes necessary; for now I'm just concerned with getting the simple meandering tests to work so I can get on with testing the new code.)
Nevertheless, a re-run of the previous test_P4_mb1_nb1_1 inputs produces the same output as before:
So...
The error above concerns adding a node with zero area. The failed assertion in the code is in ComputeVoronoiArea(), which doesn't tell me very much (although I'm wondering why the debug conditional isn't printing an error before the assertion fails?):
varea = area;
if( varea<=0.0 ) { // debug
std::cout << "Error: zero or negative varea = " << varea << std::endl;
std::cout << "Node: " << id << " " << x << " " << y << " "
<< BoundName(boundary) << std::endl;
getEdg()->TellCoords();
}
assert( varea>0.0 );
varea_rcp = 1.0/varea;
I'm more concerned about why the other interior nodes are (all) being removed, especially when the stream isn't doing much meandering. I thought maybe it had to do with the hydraulic (channel) geometry, specifically the stream widths and depths, since I am not setting them explicitly (although I'm using a channel type of 1, which uses the various HYDR_ coeffs and exponents from the input file). In looking at these, I decided it might be best to avoid bankfull events, and use the actual discharge at each node. As currently compiled, setting BANKFULLEVENT = 0 creates a fatal error:
BANKFULLEVENT: precipitation rate of a bankfull event, in m/year 0 #31.6
... ADDING INLET by resetting an existing node: Innode boundary flag =0-Non Input error: BANKFULLEVENT must be greater than zero That was a fatal error, my friend! (Set "CHILD_ABORT_ON_ERROR" to generate a crash.)
I've edited the code in tStreamNet to allow this:
// Read hydraulic geometry parameters: first, those used (potentially)
// by both "regime" and "parker" hydraulic geometry models
bankfullevent = infile.ReadItem( bankfullevent, "BANKFULLEVENT" );
bankfullevent /= SECPERYEAR;
// *SC*, 3/14/08: allow BANKFULLEVENT = 0 so that actual discharge is used in FindChanGeom()
if(0 && bankfullevent<=0.0 )
ReportFatalError(
"Input error: BANKFULLEVENT must be greater than zero" );
And in FindChanGeom() to make use of this:
// *SC*, 3/14/08: check for BANKFULLEVENT == 0 so that actual discharge is used
//qbf = cn->getDrArea()*bankfullevent;
//if( !qbf ) qbf = cn->getQ()/SECPERYEAR; // q is in m^3/s
if (bankfullevent > 0) qbf = cn->getDrArea()*bankfullevent;
else qbf = cn->getQ()/SECPERYEAR; // q is in m^3/s
width = kwds * pow(qbf, ewds);
depth = kdds * pow(qbf, edds);
rough = knds * pow(qbf, ends);
lambda = klambda * pow(qbf, elambda);
Now when I run Child, I get a different kind of error:
... tMesh::InsertNode() IC pt 3399 added at 484.975,1.71507 tMesh::InsertNode() IC pt 3400 added at 484.975,1.35727 tMesh::InsertNode() IC pt 3401 added at 484.975,0.999464 tMesh::InsertNode() IC pt 3402 added at 484.975,0.641661 tMesh::InsertNode() IC pt 3403 added at 484.975,0.283858 Mesh error: node 294 going round and round Bailing out of FlowDirs() That was a fatal error, my friend! (Set "CHILD_ABORT_ON_ERROR" to generate a crash.)
So, anyway, that's where I am at present, trying to figure out why my 'control' case isn't working right. I hope I don't have to muck around too much with the dynamic remeshing code. Once this seems to be OK, I'll run additional meandering tests comparing the (Δx, Δy) results I get from the C++ and Python code before going on to a new meandering model. Hopefully I can get this done this weekend.