JustPIC
JustPIC.jl stores material properties on particles. Once JustPIC is loaded, the geometry routines of GMG accept JustPIC particles in place of a grid. Phases and temperatures are then evaluated at the position of every particle. Interfaces are therefore resolved at particle resolution, rather than at the resolution of the numerical grid, and do not show a staircase pattern.
using GeophysicalModelGenerator, JustPIC
particles = JustPIC.init_particles(backend, nxcell, max_xcell, min_xcell, grid_vel...)
pPhases, pT = JustPIC.init_cell_arrays(particles, Val(2))
pg = particle_grid(particles; length_unit = 1m) # particle coordinates are in m
add_box!(pPhases, pT, pg; xlim = (0, 500), zlim = (-80, 0),
phase = LithosphericPhases(Layers = [20, 80], Phases = [1, 2]), T = HalfspaceCoolingTemp(Age = 80))
add_sphere!(pPhases, nothing, pg; cen = (250, 0, -100), radius = 20, phase = ConstantPhase(3))Particle grid
particle_grid collects the active particles once. Every later call reuses it and modifies the particle arrays in place, so build it after init_particles and before the particles move.
Units:
length_unitis the length of one unit of the particle coordinates, for example1m, or the characteristic length of a nondimensional model. Geometry parameters are given in km, as everywhere else inGMG.2D axes: in 2D, the second particle coordinate is the vertical one and corresponds to
zinGMG.Temperature: pass
nothinginstead of a temperature array to set only phases. Passing a thermal structureTwithout a temperature array is an error.
The following routines are supported: add_box!, add_sphere!, add_ellipsoid!, add_cylinder!, add_layer!, add_polygon!, add_plate!, add_slab!, add_volcano!, add_fault!, add_stripes!, above_surface and below_surface. The last two return a Bool array shaped like the particle CellArrays, which is false for empty particle slots.
Free surface
A GMG topography sets the free surface of a JustPIC MarkerChain (2D) or MarkerSurface (3D):
fill_chain_from_vertices!(chain, topo_line; length_unit = 1m) # topo_line: CartData of size (nx, 1, 1), e.g. from cross_section
set_topo_from_array!(surf, topo; length_unit = 1m) # topo: CartData of size (nx, ny, 1)Together with above_surface(pg, topo) for the air phase, this keeps the particles and the free surface consistent with the same topography.