PlotExt.jl cannot plot cross section

[wenrongcao]

I tried to use the PlotExt.jl module to plot a cross section, it seems the code cannot recognize field Grid in the CartData read from a timestep.

using LaMEM
using GeophysicalModelGenerator

model_output= "/Volumes/TOSHIBA EXT/LaMEM/_sb07_series/sb07_80_pl_35_50_lt/output"

timesteps,Filenames,Time = read_LaMEM_simulation(model_output)
data_cart, time = read_LaMEM_timestep(model_output,100)

The output of the CartData is :

(CartData 
    size    : (641, 2, 241)
    x       ϵ [ -2000.0 : 2000.0]
    y       ϵ [ -1.0 : 1.0]
    z       ϵ [ -1200.0 : 20.0]
    fields  : (:phase, :density, :visc_total, :visc_creep, :velocity, :pressure, :temperature, :j2_dev_stress, :j2_strain_rate, :plast_strain, :plast_dissip)
, [0.07471801])

Then I do:

include("/Users/wenrongcao/Documents/julia_codes/PlotsExt.jl")
using .PlotsExt
plot_cross_section(data_cart, field=:phase)

Error message:

type CartData has no field Grid

Stacktrace:
 [1] getproperty
   @ ./Base.jl:37 [inlined]
 [2] plot_cross_section(::CartData; field::Symbol, timestep::Nothing, dim::Int64, x::Nothing, y::Nothing, z::Nothing, aspect_ratio::Symbol)
   @ Main.PlotsExt ~/Documents/julia_codes/PlotsExt.jl:32
 [3] top-level scope

I looked at the PlotExt.jl but cannot figure out a way around to allow the code to read the proper format of x, y, z....

[boriskaus]

The way you use PlotsExt is to say:

julia> using LaMEM, Plots

does that solve the issue?

[wenrongcao]

I tried to add Plots.ji at the begining:

using Plots, LaMEM, GeophysicalModelGenerator

model_output= "/Volumes/TOSHIBA EXT/LaMEM/_sb07_series/sb07_80_pl_35_50_lt/output"

# read all avaiable time setps
timesteps,Filenames,Time = read_LaMEM_simulation(model_output)

data_cart, time = read_LaMEM_timestep(model_output,100)

include("/Users/wenrongcao/Documents/julia_codes/PlotsExt.jl")
using .PlotsExt

plot_cross_section(data_cart, field=:phase)

It gives the same error:

Adding Plots.jl plotting extensions for LaMEM
WARNING: replacing module PlotsExt.
type CartData has no field Grid

Stacktrace:
 [1] getproperty
   @ ./Base.jl:37 [inlined]
 [2] plot_cross_section(::CartData; field::Symbol, timestep::Nothing, dim::Int64, x::Nothing, y::Nothing, z::Nothing, aspect_ratio::Symbol)
   @ Main.PlotsExt ~/Documents/julia_codes/PlotsExt.jl:32
 [3] top-level scope
   @ /Volumes/TOSHIBA EXT/LaMEM/_sb07_series/sb07_80_pl_35_50_lt/jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_W0sZmlsZQ==.jl:13

[boriskaus]

The plots extension is automatically loaded if you load both LaMEM and Plots. No need to include this separately- the tutorials give done examples of that

[wenrongcao]

Anyway, the following code works without using PlotExt.jl.
I will keep checking PlotExt.jl to see what causes the issue (something to do with the access to correct format of grid data?).

using Plots, LaMEM, GeophysicalModelGenerator
model_output= "/Volumes/TOSHIBA EXT/LaMEM/_sb07_series/sb07_80_pl_35_50_lt/output"
data_cart, time = read_LaMEM_timestep(model_output,100)

x = data_cart.x.val[:,1,1]
z = data_cart.z.val[1,1,:]
phase = (data_cart.fields.phase[:,1,:])'

plot() # Create the heatmap
heatmap(x, z, phase,
        xlabel="X Axis",  # Label for the x-axis
        ylabel="Z Axis",  # Label for the z-axis
        title="Phase Map",  # Title of the plot
        color=:viridis)  # Color scheme

[boriskaus]

I think I see now what the issue is.

We are using so-called extensions in Julia which loads additional functions once the Plots package is loaded alongside LaMEM. So the function plot_cross_section is only available if you load both.

Below a MWE of this (using Julia 1.10):

julia> using LaMEM, GeophysicalModelGenerator, Plots
julia> model  = Model(Grid(nel=(16,16,16), x=[-1,1], y=[-1,1], z=[-1,1]));
julia> matrix = Phase(ID=0,Name="matrix",eta=1e20,rho=3000);
julia> sphere = Phase(ID=1,Name="sphere",eta=1e23,rho=3200);
julia> add_phase!(model, sphere, matrix)
julia> add_sphere!(model,cen=(0.0,0.0,0.0), radius=0.5);
julia> run_lamem(model,1);

Once the simulation is finished, you can plot a cross-section with:

julia> plot_cross_section(model, y=0, timestep=2, field=:phase)

Important here is that model is a model setup generated with LaMEM.jl.
It does not currently work if model is a CartData structure, even when the help may have suggested otherwise.

[boriskaus]

@wenrongcao directly plotting of CartData structs is now possible (using PR #69).
Can you please try this locally, before I release a new version of LaMEM.jl?
You can install that with:

julia>]
pkg> add LaMEM#main

[wenrongcao]

@boriskaus
I tested the given MWE and my own code on subduction, plot_cross_section in PlotExt.jl works for both the "model" and the CartData read from the output.dat (via using read_LaMEM_timestep).

PlotExt.jl is automatically included once LaMEM and Plots are used. No need to import the module again.

However, It seems the other functions in the PlotExt.jl have similar issues:

• While plot_cross_section(data_cart, y=0, field=:phase) is working, plot_cross_section_simulation(data_cart,y=0,field=:phase) does not work.

The error message is :
MethodError: no method matching plot_cross_section_simulation(::CartData; y::Int64, field::Symbol)

• My subduction model also includes topography, I did:

model_output= "/Volumes/TOSHIBA EXT/LaMEM/_sb07_series/sb07_80_pl_35_50/output"
topo_data, time = read_LaMEM_timestep(model_output,3600,surf=true)

The CartData is :

(CartData 
    size    : (641, 2, 1)
    x       ϵ [ -2000.0 : 2000.0]
    y       ϵ [ -1.0 : 1.0]
    z       ϵ [ -6.4089246 : 10.993248]
    fields  : (:topography,)
, [14.89337])

After running plot_topo(topo_data), it gives an error:

type NamedTuple has no field Topography. It seems the PlotExt.jl needs a topo parameter, which is not included in the CartData above.

Instead, I am able to read the topography by using:

topo_data, time = read_LaMEM_timestep(model_output,timestep,surf=true) 
topo = topo_data.fields.topography[:,1,1]
x = topo_data.x.val[:,1,1]
plot(x, topo, title="Topography")

My thoughts: although the provided PlotExt.jl gives a simple way to make plots, it lacks some flexibility to adjust plotting parameters such as the axis range, colormap, and labels. I think additionally providing an explicit example code to use Plots.jl to make plots can be useful. Since I have done that, I am happy to write a short document on that.

[boriskaus]

thanks. Few comments:

1. plot_cross_section_simulation actually needs info about all available timesteps; I don't think we get that from a single file, which is why using model is probably the easiest way forward.
2. you are correct that the surf option was not passed onwards. PR Fix topo plot and pass surf option to plotting timesteps #70 fixes that.
3. there was a bug in plot_topo; it should use the field topography and not Topography. I've added a fix for this in the previous PR.

The plotting routines allow passage on additional arguments, so changing axes etc. should in principle still be possible. I haven't tested this, though