Vertical velocity solver appears to be broken. E.g. test con

Vertical velocity solver is broken about thetis HOT 6 CLOSED

thetisproject commented on September 28, 2024

Vertical velocity solver is broken

from thetis.

Comments (6)

stephankramer commented on September 28, 2024

Hm, so I was looking at this at well. And initially I was blaming the vertical velocity solver as well, but note that the uv_3d does not come out the same either (look at L2 error uv), whereas that should be just the result of a projection of a constant ufl expression. And indeed if I print out the values of uv_3d they come out as:

[[1.00829158 0.30248748 0.        ]
 [1.00829158 0.30248748 0.        ]
 [1.00825087 0.30247526 0.        ]
 ...
 [1.01373007 0.30411902 0.        ]
 [0.99114942 0.29734483 0.        ]
 [0.99114942 0.29734483 0.        ]]

from thetis.

stephankramer commented on September 28, 2024

So I suspect something to do with the recent supermesh merge?

from thetis.

wence- commented on September 28, 2024

I imagine the way I rewrote projection.py screwed something up?

from thetis.

tkarna commented on September 28, 2024

I see, it indeed looks like the projection is broken on a deformed mesh. Here's a MFE

from firedrake import *

mesh2d = RectangleMesh(4, 4, 15e3, 10e3)
n_layers = 4
mesh = ExtrudedMesh(mesh2d, n_layers, 1/n_layers)
mesh.coordinates.dat.data[:, 2] *= 1 + mesh.coordinates.dat.data[:, 0]

fs_v = VectorFunctionSpace(mesh, 'DG', 1)
uv_expr = as_vector([Constant(1.0), Constant(0.3), Constant(0)])
uv = Function(fs_v, name='uv').project(uv_expr)

l2_err_uv = errornorm(uv_expr, uv)/assemble(Constant(1.0, domain=mesh2d)*dx)
print('L2 error uv {:.6e}'.format(l2_err_uv))
assert l2_err_uv < 1e-15

from thetis.

wence- commented on September 28, 2024

Ufl change is merged. So hopefully this is fixed.

from thetis.

tkarna commented on September 28, 2024

Yes, all tests pass now.

from thetis.

Recommend Projects