MODFLOW-USGS / mfusg

Public repository for MODFLOW-USG
https://water.usgs.gov/ogw/mfusg/
7 stars 4 forks source link

Possible problem with bottom averaging #12

Closed langevin-usgs closed 7 years ago

langevin-usgs commented 8 years ago

I added a new problem to my usg repository. The problem is called test046_flowdivert (https://github.com/OGW-USGS/mfusg/tree/develop/test-cmp/test046_flowdivert). The NWT model is also in there. It is a 2D model and the bottom surface rises in the middle of the domain. Constant heads cause flow from left to right, and the flow is diverted around the hump in the middle of the domain. See figure below.

This problem solves fine with mfnwt, but bottom averaging must be used. The problem solves fine with mf6, with or without bottom averaging. But the problem does not solve with mfusg, with or without damping (and thus bottom averaging). It appears that in the dry cells, the linear solver puts the head below the bottom and then the bottom averaging pulls the head back up. If I shut off the bottom averaging in USG by commenting it out (Hnew(N) = HTEMP(N)_(1.0-closebot) + Bot(N)_closebot), then it solves fine.

Joe and I talked a bit, and its not clear how best to fix this. The problem has a good mass balance, even though it is not technically converged. Perhaps convergence should be evaluated differently in this case. It is also unclear why you are using HTEMP in the bottom averaging part. It seems like that should be hnew, although that doesn't help for this problem.

sorab commented 8 years ago

Let us look at flow between 2 nodes. If head at the upstream node is below the bottom of the sill between the two nodes then the flow from the upstream node should be zero. Thus, no water goes AWAY from the node so head should not drop at that node. So, why would heads consistently drop at each iteration when head was below the bottom (the adjustment was just bringing it closer to bottom)? The culprit was line IF(EKR.LT.1.0E-7) EKR = 1.0E-7 in subroutine KR_CAL which gave that slight conductance when it should not have done so, thus sending water away from the dry node and reducing the head, which we promptly adjust back up to the bottom (almost) and that process repeats itself. Commenting out this line does the trick.

Regarding Hnew(N) = HTEMP(N)(1.0-closebot) + Bot(N)closebot), this is what was intended. If HNEW is below bottom then we are already below bottom. That is why the adjustment is to use HTEMP (which was above bottom) multiplied by a factor of 0.1 with 0.9*Bot to bring it even closer to the bottom (with reduced Kr and thus hopefully a less head reduction due to outflow from the upstream node when it is getting dry. If you use HNEW, it is already below the bottom and that is not what is intended. If head is too much below the bottom, then when there is rewetting of that node, you get a big jump in the head value (up to the bottom and then the additional physical part) and that is why I try to have it as close to the bottom as possible when CONSTANTCV is not used.

What we need in USG is what Rich has put into NWT which is to look at the bottom of the domain (all layers) and not let heads go below that whether CONSTANTCV or not.

Why it does converge when you allow heads to go below the bottom is that the numbers (Kr times the absolute K value) are small enough that the adjustment is less than the tolerance. I know you put the tolerance very small here to check this issue so I am not sure what to say here. I will play with this some more to see if I can have an answer to this.

We can discuss all this face-to-face if it is a bit confusing as there are so many little things involved here.

rniswon commented 8 years ago

Hi Sorab, Once again it is the tiniest subtleties that make such a huge difference. I now wonder why mf6 does not need any bottom correction but NWT does. I will look into it. This is why it is so good to work with you, Chris, and Joe on this stuff.

As it turns out, I will be in DC March 9-11 to present some NSF work. So maybe this will work to all get together. I imagine this is too long to wait to get to the bottom of this specific issue but it may be a good opportunity none the less.

Rich

On Dec 8, 2015, at 5:25 PM, Sorab notifications@github.com wrote:

Let us look at flow between 2 nodes. If head at the upstream node is below the bottom of the sill between the two nodes then the flow from the upstream node should be zero. Thus, no water goes AWAY from the node so head should not drop at that node. So, why would heads consistently drop at each iteration when head was below the bottom (the adjustment was just bringing it closer to bottom)? The culprit was line IF(EKR.LT.1.0E-7) EKR = 1.0E-7 in subroutine KR_CAL which gave that slight conductance when it should not have done so, thus sending water away from the dry node and reducing the head, which we promptly adjust back up to the bottom (almost) and that process repeats itself. Commenting out this line does the trick.

Regarding Hnew(N) = HTEMP(N)(1.0-closebot) + Bot(N)closebot), this is what was intended. If HNEW is below bottom then we are already below bottom. That is why the adjustment is to use HTEMP (which was above bottom) multiplied by a factor of 0.1 with 0.9*Bot to bring it even closer to the bottom (with reduced Kr and thus hopefully a less head reduction due to outflow from the upstream node when it is getting dry. If you use HNEW, it is already below the bottom and that is not what is intended. If head is too much below the bottom, then when there is rewetting of that node, you get a big jump in the head value (up to the bottom and then the additional physical part) and that is why I try to have it as close to the bottom as possible when CONSTANTCV is not used.

What we need in USG is what Rich has put into NWT which is to look at the bottom of the domain (all layers) and not let heads go below that whether CONSTANTCV or not.

Why it does converge when you allow heads to go below the bottom is that the numbers (Kr times the absolute K value) are small enough that the adjustment is less than the tolerance. I know you put the tolerance very small here to check this issue so I am not sure what to say here. I will play with this some more to see if I can have an answer to this.

We can discuss all this face-to-face if it is a bit confusing as there are so many little things involved here.

— Reply to this email directly or view it on GitHub https://github.com/OGW-USGS/mfusg/issues/12#issuecomment-163077044.

sorab commented 8 years ago

One more subtlety (when using CONSTANTCV) and a single layer and starting with heads below the bottom and if there is recharge (or injection well) then without that small number the whole row of the matrix becomes zero (but recharge is on RHS). We take care of that in GLO2SMS1AP before entering the solvers. We currently do the following: C3c---------TAKE CARE OF ZERO ROW DIAGONAL ADIAG = ABS(AMAT(IA(N))) IF(ADIAG.LT.1.0E-15)THEN AMAT(IA(N)) = 1.0E06 RHS(N) = RHS(N) + HNEW(N)*1.0E06 ENDIF

We should update RHS(N) as per above, only if RHS(N) is zero (abs value less than 1e15). That will not shock the system for the above case.

langevin-usgs commented 8 years ago

So what is the entire block of code that you have for C3c? Can you paste it here?

Also, it seems a little strange to tie all this to CONSTANTCV, doesn't it? I don't think someone would expect the presence or absence of CONSTANTCV to have any bearing on a one-layer model.

BTW, my testing so far has showed overall good things for setting the "magic number" from 1e-7 to 0. One of the USG models no longer converges (../examples/03A_conduit_unconfined/ex3A.nam), but that can probably be fixed. In MF6, this change is also working so far, but more testing is underway. Joe had also implemented a search for each node that finds the lowest bottom elevation beneath it (although not recursively yet) and is using that for the bottom averaging. The HTEMP/HNEW issue makes sense according to your explanation.

sorab commented 8 years ago

Block of code is below: C3c---------TAKE CARE OF ZERO ROW DIAGONAL ADIAG = ABS(AMAT(IA(N))) IF(ADIAG.LT.1.0E-15)THEN AMAT(IA(N)) = 1.0E06 IF(ABS(RHS(N)) .GT. 1.0E-15) 1 RHS(N) = RHS(N) + HNEW(N)*1.0E06 ENDIF The issue is all with the matrix. If you have CONSTANTCV and multiple layers then the matrix has at least one non-zero for a dry node (the K between layers). However, with only 1 layer, the entire row can be zero. Without CONSTANTCV, the head will not be below the bottom of the cell so even if the flow terms are all zero, the Newton term will not be zero (epsilon shift will make the cell slightly saturated so the derivative is not zero). However, with CONTANTCV the head can drop below the bottom. Then, with an epsilon shift for the Newton derivative, the head will very likely be STILL below the bottom (still dry with Kr=0) so the derivative is zero and that can cause issues.

sorab commented 8 years ago

Please ignore above lines. See lines below.

C3c---------TAKE CARE OF ZERO ROW DIAGONAL ADIAG = ABS(AMAT(IA(N))) IF(ADIAG.LT.1.0E-15)THEN AMAT(IA(N)) = 1.0 IF(ABS(RHS(N)) .LT. 1.0E-15) 1 RHS(N) = HNEW(N) ENDIF

The 1 on diagonal of AMAT is random, just to kick things up. If you use a smaller number, the jump at the first iteration would be larger (and vice versa). If the jump is too small and things still stay dry then you will be taking many iterations to wet the cell.

sorab commented 8 years ago

Sorry, one more time on C3c. I am including the provision for the number to be different. C3c---------TAKE CARE OF ZERO ROW DIAGONAL ADIAG = ABS(AMAT(IA(N))) IF(ADIAG.LT.1.0E-15)THEN
ANUMBER = 1.0 AMAT(IA(N)) = ANUMBER IF(ABS(RHS(N)) .LT. 1.0E-15) 1 RHS(N) = HNEW(N) * ANUMBER ENDIF

This does not happen too often (I have seen it in only a couple of problems) but it is good to have this flexibility.

sorab commented 8 years ago

I have been running tests on this. What I had previously was correct for this part. More flexibility can be got by not using 1e6 though, as per below.

C3c---------TAKE CARE OF ZERO ROW DIAGONAL ADIAG = ABS(AMAT(IA(N))) IF(ADIAG.LT.1.0E-15)THEN
ANUMBER = 1.0 AMAT(IA(N)) = ANUMBER RHS(N) = RHS(N) + HNEW(N) * ANUMBER ENDIF

Sorry for all the back and forth. I should have done all this before shooting off messages here.

langevin-usgs commented 8 years ago

For cleanliness, I suggest the following slight tweaks for double precision. The only remaining problem I see for this change is that one of the unconfined conduit problems that we distribut with mfusg no longer converges (03A_conduit_unconfined/ex3A.nam).

C3c---------TAKE CARE OF ZERO ROW DIAGONAL ADIAG = ABS(AMAT(IA(N))) IF(ADIAG.LT.1.0D-15)THEN ANUMBER = 1.D0 AMAT(IA(N)) = ANUMBER RHS(N) = RHS(N) + HNEW(N) * ANUMBER ENDIF

langevin-usgs commented 7 years ago

This was fixed in https://github.com/OGW-USGS/mfusg/commit/2f397cba708034205d983a444aac56c28ccc0517. Closing this issue.