How much effect will this change have on the multi-phase solver formulation?

Concerning multiphase:
We originally switched it for the sake of doing abl+waves. It helps avoid advecting the huge discontinuity in rho*theta. However, that doesn't mean it's right. I was mainly concerned about maintaining the temperature-based dynamics in the atmosphere and keeping things stable at the phase interface, while understanding the treatment at the interface might not be realistic. Two thoughts.

1. If it is indeed physically correct to model two-phase temperature transport by a single theta equation, then we can address the advection of rhotheta the same way we do it for rhou. This shouldn't be that hard to set up, I'll just take some code from icns_advection.H and put it in AdvOp_Godunov.H (something like that).
2. I have doubts that two-phase temperature can be modeled by a single equation. If you consider temperature from air and temperature from water coming into a single cell, the equilibrium temperature of the cell shouldn't be equal to just the mass-weighted average of the two temperatures; the specific heat of each material would have to be taken into account. This points to (probably) the fact that we'd need to be solving for a single energy equation from which we can calculate the temperature. This would be a new thing for AMR-Wind, seems like a harder approach.

The complexity here is why I say we should just revert back to multiply_rho = true. For multiphase flows, we can implement 1., though we might not be able to trust the thermal exchange at the interface. It may be a good stand-in for the full approach, since water is has both a much higher density and a much higher heat capacity than air. Also, I might be wrong about needing an energy equation; there may be a way to do it with a single temperature equation + source terms.

I checked the implementation of VOF based methods in OpenFOAM and temperature equation is not included in the calculations. The temperature inclusion makes it necessary to incorporate additional physics: https://github.com/OpenFOAM/OpenFOAM-6/tree/master/applications/solvers/multiphase/multiphaseEulerFoam. I think an easier option for us is to specify a temperature profile within the water region and apply a forcing to ensure that the temperature within the water remains horizontally homogeneous to avoid any artifacts at the interface. We also need to ensure consistency of the temperature profile with the air temperature at the interface to avoid any abrupt jumps.

I just added a bunch of stuff to make it work in #1498, see here: marchdf#2

The key concept is that the rho part of rhoT is transported the same way as rho, which makes things stable. There's a little more complexity in there too, but it's all been tested in the mass-momentum consistency context, which we've published. This will make sure things don't blow up, but I'm still doubtful that this is a physically realistic equation to solve for both phases. We'll figure that out later. One thing to note: this still assumes that the density of each phase is constant. Doing multiphase + anelastic would require more coding.

I did confirm that it works for abl_multiphase_laminar. Because I had to change something in incflo_advance, we need to make sure that the other tests don't end up with new diffs as a result.