implement "hydrostatic compression" mass term #2012
Conversation
tjhei
force-pushed
the
hydrostatic_compression_formulation
branch
2 times, most recently
from
af434cd
to
2e8ddf0
Compare
| @@ -120,7 +120,7 @@ namespace aspect | |||
| * includes this term implicitly in the matrix, | |||
| * which is therefore not longer symmetric. | |||
| * This class approximates this term as | |||
| * $ - \nabla \mathbf{u} - \frac{1}{\rho^{\ast}} * \frac{\partial rho{^\ast}}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u} = 0$ | |||
| * $ - \nabla \cdot \mathbf{u} - \frac{1}{\rho^{\ast}} * \frac{\partial rho{^\ast}}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u} = 0$ | |||
There was a problem hiding this comment.
The cdot is also missing in the documentation of the StokesReferenceDensityCompressibilityTerm in line 103.
|
|
||
| /** | ||
| * This class assembles the right-hand-side term of the Stokes equation | ||
| * that is caused by the compressibility in the mass conservation equation. |
There was a problem hiding this comment.
mention that it includes both temperature and pressure? In the old pull request, we had
caused by the compression based on the changes in pressure and temperature in the mass conservation equation
| @@ -447,7 +447,11 @@ namespace aspect | |||
| internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>& > (data_base); | |||
|
|
|||
| // assemble RHS of: | |||
| // - div u = 1/rho * drho/dp rho * g * u | |||
| // - div \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial p} \rho \mathbf{g} \cdot \mathbf{u} | |||
There was a problem hiding this comment.
The * shouldn't be there.
|
|
||
| // assemble RHS of: | ||
| // - div \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial p} \rho \mathbf{g} \cdot \mathbf{u} | ||
| // + \frac{1}{\rho} * \frac{\partial rho}{\partial T} \nabla T \cdot \mathbf{u} |
There was a problem hiding this comment.
And the alpha term has the wrong sign (it should be - \frac{1}{\rho} \frac{\partial rho}{\partial T} \nabla T \cdot \mathbf{u}). See the implementation (which has the correct sign).
tests/hydrostatic_compression.cc
Outdated
| "error calculation", | ||
| "A postprocessor that compares the numerical solution to the analytical " | ||
| "solution derived for compressible mantle convection in a 2D box as described " | ||
| "in the manuscript and reports the error.") |
There was a problem hiding this comment.
Which manuscript? This test seems to have some kind of analytical solution, so it would be nice to have a bit more documentation about that.
| RMS, max velocity: 29.3 m/s, 42.2 m/s | ||
| Pressure min/avg/max: -0.002532 Pa, 0.4959 Pa, 0.989 Pa | ||
| Max and min velocity along boundary parts: 42.03 m/s, 38.83 m/s, 42.03 m/s, 38.83 m/s, 42.03 m/s, 6.236 m/s, 42.01 m/s, 6.067 m/s | ||
| Errors ndofs= 3556 u_L2= 1.001381e-02 p_L2= 1.151619e-02 |
There was a problem hiding this comment.
Should we look if this converges?
There was a problem hiding this comment.
it does and the results are in the paper manuscript
tjhei
force-pushed
the
hydrostatic_compression_formulation
branch
from
2e8ddf0
to
e39f890
Compare
|
okay, ready for rereview |
|
|
||
| // Compared to the manual, this term seems to have the wrong sign, but | ||
| // this is because we negate the entire equation to make sure we get | ||
| // -div(u) as the adjoint operator of grad(p) |
There was a problem hiding this comment.
Please remove the duplicate documentation (copy and paste problem). Note that the first version is missing the word "compared".
tjhei
force-pushed
the
hydrostatic_compression_formulation
branch
from
e39f890
to
d3ee16b
Compare
|
👍 |
|
Nice! I suppose that the testcase has an analytic solution you double checked? |
Yes, I see optimal (cubic) convergence in |u|_0 for the test specifically constructed to have a hydrostatic density. We wanted to have this PR merged to continue our experiments. I think we are planning to do more experiments before we think about advertising this inside ASPECT. I assume this should become the new default formulation but it is too early for now. |
supersedes #1407