Nondimensionalization of Geothermal Gradient

[aelligp]

I've ran into a nondimesionalization problem during testing of my code.

I want to calculate the geothermal gradient on a grid that has a topography. This means that the 0.0km mark should start at topography. After doing so, I wanted to add the gradient in non dimensional numbers.

However the non dimensional numbers do not work correctly when dimensionalized again. See screenshot.

Is there a way to add a non dimensional term for K/km or C/km?

[boriskaus]

Can you add a full MWE (minimal working example)? Screenshots are only partially helpful.

[aelligp]

using GeoParams

CharDim     = GEO_units(length=40km, viscosity=1e20Pa*s);

Depth       = Array(0km:1km:10km);
Depth_nondim= nondimensionalize(Depth,CharDim);
Geotherm    = nondimensionalize(30K, CharDim)/nondimensionalize(1km, CharDim)
Geotherm_C  = nondimensionalize(30C, CharDim)/nondimensionalize(1km, CharDim)

Gradient    = nondimensionalize(273.15K,CharDim) .+ Geotherm * Depth_nondim;
Temp_K_dim  = dimensionalize(Gradient, K, CharDim)

Gradient_C  = nondimensionalize(0C,CharDim) .+ Geotherm_C * Depth_nondim;
Temp_C_dim  = dimensionalize(Gradient_C, C, CharDim)

and these are the results from this:

11-element Vector{Quantity{Float64, 𝚯, Unitful.FreeUnits{(K,), 𝚯, nothing}}}:
             273.15 K
             303.15 K
             333.15 K
             363.15 K
             393.15 K
  423.1499999999999 K
             453.15 K
 483.15000000000003 K
             513.15 K
             543.15 K
             573.15 K


11-element Vector{Quantity{Float64, 𝚯, Unitful.FreeUnits{(K,), 𝚯, Unitful.Affine{-5463//20}}}}:
                0.0 °C
             303.15 °C
              606.3 °C
  909.4499999999999 °C
             1212.6 °C
 1515.7499999999995 °C
 1818.8999999999992 °C
 2122.0499999999993 °C
 2425.1999999999994 °C
            2728.35 °C
 3031.4999999999995 °C

[boriskaus]

ok, I see. Here a more compact version of the same:

julia> using GeoParams
julia> CharDim    = GEO_units(length=40km, viscosity=1e20Pa*s);
julia> GeoTherm = GeoUnit(30K/1km)
GeoUnit{dimensional, K km⁻¹·⁰}, 
30.0
julia> GeoTherm_nd = nondimensionalize(GeoTherm, CharDim)
GeoUnit{nondimensional, K km⁻¹·⁰}, 
0.942544083572242

julia> GeoTherm_dim = dimensionalize(GeoTherm_nd, CharDim)
julia> @assert GeoTherm_dim == GeoTherm.  # all is fine

julia> length = GeoUnit(1km)
julia> length_nd = nondimensionalize(length, CharDim)
GeoUnit{nondimensional, km}, 
0.025

julia> T = length*GeoTherm.  # correct answer
GeoUnit{dimensional, K}, 
30.0

julia> T_nd = length_nd*GeoTherm_nd. # incorrect units are stored (should be {nondimensional, K})
GeoUnit{nondimensional, }, 
0.02356360208930605

So there is a problem when we multiply two non-dimensional GeoUnit parameters.

[albert-de-montserrat]

As far as I can tell from here, units are thrown away.

[boriskaus]

As far as I can tell from here, units are thrown away.

yup, that's no good.

[albert-de-montserrat]

This may work?

function foo(x::GeoUnit{T1,U1}, y::GeoUnit{T2,U2}) where {T1,T2,U1,U2}
    return GeoUnit(*(UnitValue(x), UnitValue(y)) * (x.unit*y.unit))
end

In [21]: foo(length_nd, GeoTherm_nd)
GeoUnit{dimensional, K},
0.02356360208930605

[boriskaus]

yes that should do it; now it needs to be generalized for all operations. Addition/subtraction can only be done if they have the same time; multiplication/division should work for different units as well

[albert-de-montserrat]

I missed it, but it returns GeoUnit{dimensional, K} instead of nondimensional

[boriskaus]

yes missed that too; should be nondimensional. ideally there should be a check to ensure that the types are the same (either dimensional or nondimensional)