search

nextsimhub / nextsimdg

neXtSIM_DG : next generation sea-ice model with DG
https://nextsim-dg.readthedocs.io/en/latest/?badge=latest
Apache License 2.0
10 stars 13 forks source link

Dynamics with spherical coordinates #548

Closed timspainNERSC closed 1 month ago

timspainNERSC commented 1 month ago

Dynamics with spherical coordinates

Fixes #463

Change Description

Clean up the coordinate handling such that spherical coordinates are passed to the dynamics.

Clean up the scaling within the dynamics such that running the 25 km Arctic Ocean regional model works with both Cartesian coordinates and spherical coordinates.

Test Description

Manual integration testing: the 25 km Arctic Ocean regional model runs with physically reasonable results for both mEVP and BBM.

Documentation Impact

The model now works as documented.

timspainNERSC commented 1 month ago

Many of the changes are the result of running clang-format on files in the dynamics for the first time.

winzerle commented 1 month ago

Hi,

In spherical coordinates, the gradient must be scaled by 1/R as all the numerics is set up on a reference element of size 1 and then scaled to the real one. But it looks like this is already correctly done in ParametricMap when the fields divS1, divS2 and divM are assembled. So the stressScale should indeed be removed.

This is what they mean:

divS1 is the

d_x Phi * Psi

Where Phi refers to the CG-part and Psi to the DG part. They implement

V = div (S)

But written as

(V, Phi) = (div S, Phi) = - (S, d_x Phi)

divS2 is the same for y.

In spherical coordinates the divergence is scaled by an 1/R. This however is already done in ParametricMap. Also for divM, which is the extra term arising in spherical coordinates coming from the derivative of the unit vectors.

I checked my spherical-coordinates benchmark and it looks good, and the shear is much more reasonable than in the case including the wrong stress-scaling.

Somewhere I put down all the stuff regarding the background higher order implementation of spherical coordinates which was rather tricky. As soon as I find the time, I’ll update it and we can put it as technical explanation somewhere in the repo.

Thomas

On 13. May 2024, at 10:11, Einar Örn Ólason @.***> wrote:

@einola commented on this pull request. Let's talk about the stressScale stuff today.In core/src/modules/DynamicsModule/MEVPDynamics.cpp:

@@ -16,6 +16,9 @@

namespace Nextsim {

+// Degrees to radians as a hex float

Why is this not needed in BBMDynamics.cpp? In dynamics/src/include/VPCGDynamicsKernel.hpp:

s11 = stress[I11]; s12 = stress[I12]; s22 = stress[I22];

  • double stressScale = 1.0; // 2nd-order Stress term has different scaling with the EarthRadius
  • if (smesh->CoordinateSystem == Nextsim::SPHERICAL)
  • stressScale = 1.0 / Nextsim::EarthRadius / Nextsim::EarthRadius;
  • double stressScale = 1.0; // Stress scale of 1 correctly scales stress in both CARTESIAN

That's ... interesting? Did you ask Thomas why there's this stressScale to begin with? If this is really the case, then we should just get rid of the stressScale variable. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because your review was requested.Message ID: @.***>

einola commented 1 month ago

Thank you for the explanation, Thomas. That makes sense :)