Read timesteps back into LaMEM

If you want to quantitatively do something with the results, there is an easy way to read the output of a LaMEM timestep back into julia. All routines related to that are part of the LaMEM.IO module.

julia> using LaMEM

You can first read the *.pvd file in the directory to see which timesteps are available:

julia> FileName="FB_multigrid"
julia> DirName ="test"
julia> Timestep, Filenames, Time = read_LaMEM_simulation(FileName, DirName)
([0, 1], ["Timestep_00000000_0.00000000e+00/FB_multigrid.pvtr", "Timestep_00000001_6.72970343e+00/FB_multigrid.pvtr"], [0.0, 6.729703])

We can read a particular timestep (say 1) with:

julia> data, time = read_LaMEM_timestep(FileName, 1, DirName)
(CartData 
    size    : (33, 33, 33)
    x       ϵ [ 0.0 : 1.0]
    y       ϵ [ 0.0 : 1.0]
    z       ϵ [ 0.0 : 1.0]
    fields  : (:phase, :visc_total, :visc_creep, :velocity, :pressure, :strain_rate, :j2_dev_stress, :j2_strain_rate)
  attributes: ["note"]
, [6.729703])

The output is in a CartData structure (as defined in GeophysicalModelGenerator).

If you do not indicate a directory name (DirName) it'll look in your current directory. The default above will load the main LaMEM simulation output. Alternatively, you can also load the phase information by specify the optional keyword phase=true:

julia> data, time = read_LaMEM_timestep(FileName, 1, DirName, phase=true)
(CartData 
    size    : (96, 96, 96)
    x       ϵ [ 0.0052083334885537624 : 0.9947916269302368]
    y       ϵ [ 0.0052083334885537624 : 0.9947916269302368]
    z       ϵ [ 0.0052083334885537624 : 0.9947916269302368]
    fields  : (:phase,)
  attributes: ["note"]
, [6.729703])

In the same way, you can load the internal free surface with surf=true (if that was saved), or passive tracers (passive_tracers=true). If you don't want to load all the fields in the file back to julia, you can check which fields are available:

julia> read_LaMEM_fieldnames(FileName, DirName)
("phase [ ]", "visc_total [ ]", "visc_creep [ ]", "velocity [ ]", "pressure [ ]", "strain_rate [ ]", "j2_dev_stress [ ]", "j2_strain_rate [ ]")

and load only part of those:

julia> data, time = read_LaMEM_timestep(FileName, 1, DirName, fields=("phase [ ]", "visc_total [ ]","velocity [ ]"))
(CartData 
    size    : (33, 33, 33)
    x       ϵ [ 0.0 : 1.0]
    y       ϵ [ 0.0 : 1.0]
    z       ϵ [ 0.0 : 1.0]
    fields  : (:phase, :visc_total, :velocity)
  attributes: ["note"]
, [6.729703])

Removing the output of a run

The results of an earlier run can be cleared with clean_directory, which takes either the model or the name of a directory:

julia> clean_directory(model)              # the output directory of the model
julia> clean_directory("my_output_dir")    # or a directory by name

This removes the Timestep_* directories and the *.pvd/*.vts files, but keeps the directory itself along with the *.dat input file, the markers and any logfile, so the model can be run again. Use rm(model.Output.out_dir, recursive=true) if you want the directory gone as well.

Clearing the output is worth doing when you change a setup and then look at it with view_model: without it, the viewer finds the timesteps of the previous run and offers their fields rather than those of the new setup.