Temperature structure of LinearTemp()

[aelligp]

When adding e.g. a cylinder to a grid with a volcano and applying a linear temperature distribution, I dont get the correct Ttop and Tbot. LinearTemp seems to assume the Ttop to be at 0.

MWE:

    Lx = Ly = 40
    nx = ny = nz = 32
    x = range(-Lx / 2, Lx / 2, nx);
    y = range(-Ly / 2, Ly / 2, ny);
    z = range(-30, 10, nz);
    Grid = CartData(xyz_grid(x,y,z));

    Phases = fill(5, nx, ny, nz);
    Temp = fill(20.0, nx, ny, nz);

    lith = LithosphericPhases(Layers=[15 45 100], Phases=[1 2 3])

    add_box!(Phases, Temp, Grid;
        xlim=(-100, 100),
        ylim=(-400, 400.0),
        zlim=(-110.0, 0.0),
        phase = lith,
        T = HalfspaceCoolingTemp(Age=20)
    )

    add_volcano!(Phases, Temp, Grid;
        volcanic_phase  = 1,
        center          = (0, 0, 0),
        height          = 4,
        radius          = 5,
        crater          = 0.5,
        base            = 0.0,
        background      = nothing,
        T               = HalfspaceCoolingTemp(Age=20)
    )

    add_ellipsoid!(Phases, Temp, Grid;
        cen    = (0, 0, -5),
        axes   = (2.5, 2.5, 2.5/2),
        phase  = ConstantPhase(3),
        T      = ConstantTemp(T=1000)
    )
    add_cylinder!(Phases, Temp, Grid;
        base = (0, 0, -3.75),
        cap  = (0, 0, 4.0),
        radius = 1.0,
        phase  = ConstantPhase(3),
        T      = LinearTemp(Ttop=20, Tbot=1000)
    )

   Grid = addfield(Grid,(; Phases, Temp))
    write_paraview(Grid, "MWE")

Here a high resolution output of it:

[aelligp]

The same thing happens with an ellipsoid below the surface.