Computation of gradients in a compressible SingleStack

[yairchn]

The EDMF development has been struggling with the computation of gradients (in the z direction) in the 1D model.
To isolate the issues I have created an example branch 'yc/gradient_exmaple':

https://github.com/CliMA/ClimateMachine.jl/tree/yc/gradient_exmaple

In this case I am running a compressible singlestack for 100sec in a Stable BL and plotting both the u velocity component and its vertical gradient:

[yairchn]

This can be reproduced by:

1. running
   using Revise; include("test/Atmos/EDMF/stable_bl_coupled_edmf_an1d.jl")

2. after the simulation is done pasting in Julia RPEL
   """"
   const clima_dir = dirname(dirname(pathof(ClimateMachine)));
   include("test/Atmos/EDMF/compute_mse.jl")
   include(joinpath(clima_dir, "docs", "plothelpers.jl"));
   dons_arr = get_dons_arr(diag_arr)
   output_dir = joinpath(clima_dir, "output")
   z = Array(get_z(solver_config.dg.grid; rm_dupes = true))

export_plot(
z,
time_data,
dons_arr,
("∇flux_turbconv_∇u_3",),
joinpath(output_dir, "dudz.png");
xlabel = "∇u (1 / s)",
ylabel = "z (m)",
time_units = "(seconds)",
lw = 2,
palette = :tab10,
)

export_plot(
z,
time_data,
dons_arr,
("prog_ρu_1",),
joinpath(output_dir, "ρu.png");
xlabel = "ρu (kg / m^2 s)",
ylabel = "z (m)",
time_units = "(seconds)",
lw = 2,
palette = :tab10,
)
""""

[yairchn]

The behavior of 'dudz' in the plot above for 20<z<100 is hard to connected with the profile of u.
Examining the 300m level it seems that the DG I trying to keep the derivative smooth (continuous) but it does not produces a desirable outcome.

[jkozdon]

I think that the derivative is correct. The fields are the problem. I did a quick and dirty finite difference of the field u and got similar results.

for i = 1:length(dons_arr)
  dons_arr[i]["prog_u"] = dons_arr[i]["prog_ρu_1"] ./ dons_arr[i]["prog_ρ"]
  Δz = diff(dons_arr[i]["z"])
  Δu = diff(dons_arr[i]["prog_u"])
  Δρu = diff(dons_arr[i]["prog_ρu_1"])
  Δρ = diff(dons_arr[i]["prog_ρ"])
  Δu_Δz = Δu ./ Δz
  Δρu_Δz = Δρu ./ Δz
  Δρ_Δz = Δρ ./ Δz
  dons_arr[i]["Δu_Δz"] = [Δu_Δz[1]; (Δu_Δz[1:end-1] + Δu_Δz[2:end]) / 2; Δu_Δz[end]]
  dons_arr[i]["Δρu_Δz"] = [Δρu_Δz[1]; (Δρu_Δz[1:end-1] + Δρu_Δz[2:end]) / 2; Δρu_Δz[end]]
  dons_arr[i]["Δρ_Δz"] = [Δρ_Δz[1]; (Δρ_Δz[1:end-1] + Δρ_Δz[2:end]) / 2; Δρ_Δz[end]]
end
export_plot(
            z,
            time_data,
            dons_arr,
            ("Δu_Δz",),
            joinpath(output_dir, "Δu_Δz.png");
            xlabel = "Δu_Δz",
            ylabel = "z (m)",
            time_units = "(seconds)",
            lw = 2,
            palette = :tab10,
           )
export_plot(
            z,
            time_data,
            dons_arr,
            ("Δρu_Δz",),
            joinpath(output_dir, "Δρu_Δz.png");
            xlabel = "Δρu_Δz",
            ylabel = "z (m)",
            time_units = "(seconds)",
            lw = 2,
            palette = :tab10,
           )
export_plot(
            z,
            time_data,
            dons_arr,
            ("Δρ_Δz",),
            joinpath(output_dir, "Δρ_Δz.png");
            xlabel = "Δρ_Δz",
            ylabel = "z (m)",
            time_units = "(seconds)",
            lw = 2,
            palette = :tab10,
           )

which produces the following plots

[yairchn]

@jkozdon thanks so much. That really helps me isolating where the problem is.
I will work on it from here and close this issue

[yairchn]

I am re-opening the issue with a better view of the problem following @jkozdon's comment.
I can reproduce this issue with the ConstantKinematicViscosity(0.1) model with no EDMF (using 'NoTurbConv') flag.
I pushed these changes so that running the case above on this branch will now run the above settings and pasting this

""""""""""""""
const clima_dir = dirname(dirname(pathof(ClimateMachine)));
include("test/Atmos/EDMF/compute_mse.jl")
include(joinpath(clima_dir, "docs", "plothelpers.jl"));
dons_arr = get_dons_arr(diag_arr)
output_dir = joinpath(clima_dir, "output")
z = Array(get_z(solver_config.dg.grid; rm_dupes = true))

export_plot(
z,
time_data,
dons_arr,
("∇flux_u_3",),
joinpath(output_dir, "dudz_gm.png");
xlabel = "∇u (1 / s)",
ylabel = "z (m)",
time_units = "(seconds)",
lw = 2,
palette = :tab10,
)

export_plot(
z,
time_data,
dons_arr,
("prog_ρu_1",),
joinpath(output_dir, "ρu.png");
xlabel = "ρu (kg / m^2 s)",
ylabel = "z (m)",
time_units = "(seconds)",
lw = 2,
palette = :tab10,
)
"""""""""""

in the RPEL will reproduce the following plots of the vertical gradient of the atmosmodel velocities:

[yairchn]

it appears that the problem is more generally a problem with the way 'flux_second_order' (which is the only flux operating here).

[jkozdon]

I'm not convinced that there is problem with the dg code here. My hunch is that it is something like Runge's phenomena you are seeing.

I think that the most helpful thing would be to write down the exact mathematically equations you are trying to solve including the domain and the boundary conditions. At this point I have little idea how all the driver and atmos model actually define a problem, so it's hard to make suggestions.

Once a simplified problem is written down, there are simpler DG code that you could implement the problem in to explore what you are seeing and possible fixes.