added the possibility to create a CartGrid with dimensional values & non-dimensionalise it

Author: boriskaus

src/data_types.jl

@@ -962,28 +962,45 @@ end
 
 """
 
-Creates a 1D, 2D or 3D cartesian grid of given size. Grid can be created by defining the size and either the `extent` (length) of the grid in all directions, or by defining start & end points 
+    Grid = CreateCartGrid(; size=(), x = nothing, z = nothing, y = nothing, extent = nothing, CharDim = nothing)
 
-Spacing is assumed to be constant
+Creates a 1D, 2D or 3D cartesian grid of given size. Grid can be created by defining the size and either the `extent` (length) of the grid in all directions, or by defining start & end points (`x`,`y`,`z`). 
+If you specify `CharDim` (a structure with characteristic dimensions created with `GeoParams.jl`), we will nondimensionalize the grd before creating the struct.
+
+Spacing is assumed to be constant in a given direction
+
+This can also be used for staggered grids, as we also create 1D vectors for the central points. The points you indicate in `size` are the corner points.
+
+Note: since this is mostly for solid Earth geoscience applications, the second dimension is called z (vertical)
 
-Note: since this is mostly for Solid Earth geoscience applications, the second dimension is called z (vertical)
 
 # Examples
 ====
 
+A basic case with non-dimensional units:
 ```julia
-Grid = CreateCartGrid(size=(10,20),x=(0.,10), z=(2.,10))
+julia> Grid = CreateCartGrid(size=(10,20),x=(0.,10), z=(2.,10))
 Grid{Float64, 2} 
            size: (10, 20) 
          length: (10.0, 8.0) 
          domain: x ∈ [0.0, 10.0], z ∈ [2.0, 10.0] 
  grid spacing Δ: (1.1111111111111112, 0.42105263157894735) 
 ```
+
+An example with dimensional units:
+```julia
+julia> CharDim = GEO_units()
+julia> Grid    = CreateCartGrid(size=(10,20),x=(0.0km, 10km), z=(-20km, 10km), CharDim=CharDim)
+
+```
+
+
 """
 function CreateCartGrid(;
     size=(),
      x = nothing, z = nothing, y = nothing,
-     extent = nothing
+     extent = nothing,
+     CharDim = nothing
 )
 
     if isa(size, Number)
@@ -1023,6 +1040,18 @@ function CreateCartGrid(;
     end
     Xₙ  = X₁ .+ L  
     Δ   = L ./ (N .- 1)       
+    
+    # nondimensionalize 
+    if !isnothing(CharDim)
+        X₁, Xₙ, Δ, L    = GeoUnit.(X₁), GeoUnit.(Xₙ), GeoUnit.(Δ),  GeoUnit.(L)
+        
+        X₁              = ntuple( i -> nondimensionalize(X₁[i], CharDim), dim)
+        Xₙ              = ntuple( i -> nondimensionalize(Xₙ[i], CharDim), dim)
+        Δ               = ntuple( i -> nondimensionalize(Δ[i],  CharDim), dim)
+        L               = ntuple( i -> nondimensionalize(L[i],  CharDim), dim)
+
+        X₁, Xₙ, Δ, L    = NumValue.(X₁), NumValue.(Xₙ), NumValue.(Δ), NumValue.(L)
+    end
 
     # Generate 1D coordinate arrays of vertexes in all directions
     coord1D=()
@@ -1041,6 +1070,8 @@ function CreateCartGrid(;
 
 end
 
+
+
 # view grid object
 function show(io::IO, g::CartGrid{FT, DIM}) where {FT, DIM}
 

test/test_data_types.jl

@@ -151,6 +151,10 @@ Data     =   Z
 Data_setC =   CartData(X,Y,Z, (FakeData=Data,Data2=Data.+1.))
 @test sum(abs.(Value(Data_setC.x))) ≈ 2310.0km
 
+
+
+
+
 # Convert from CartData -> UTMData
 Data_set4 = Convert2UTMzone(Data_setC, proj)
 

test/test_setup_geometry.jl

@@ -79,6 +79,8 @@ Phases              =   Compute_Phase(Phases, Temp, Grid, LP);
 
 # test AboveSurface with the Grid object
 Grid        =   CreateCartGrid(size=(10,20,30),x=(0.,10), y=(0.,10), z=(-10.,2.))
+@test Grid.Δ[2] ≈ 0.5263157894736842
+
 Temp        =   ones(Float64, Grid.N...)*1350;
 Phases      =   zeros(Int32,  Grid.N...);
 
@@ -89,3 +91,8 @@ ind         =   AboveSurface(Grid, Topo_cart);
 
 ind         =   BelowSurface(Grid, Topo_cart);
 @test sum(ind[1,1,:]) == 25
+
+# Create Grid & nondimensionalize it
+CharDim     =   GEO_units();
+Grid        =   CreateCartGrid(size=(10,20,30),x=(0.0km,10km), y=(0.0km, 10km), z=(-10.0km, 2.0km), CharDim=CharDim)
+@test Grid.Δ[2] ≈ 0.0005263157894736842