Skip to content

Commit c56c877

Browse files
boriskausboriskaus
committed
Merge remote-tracking branch 'origin/main' into main
2 parents 592d903 + 3c82495 commit c56c877

File tree

5 files changed

+375
-55
lines changed

5 files changed

+375
-55
lines changed

‎docs/src/man/lamem.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@ GeophysicalModelGenerator.AddBox!
2020
GeophysicalModelGenerator.AddSphere!
2121
GeophysicalModelGenerator.AddEllipsoid!
2222
GeophysicalModelGenerator.AddCylinder!
23+
GeophysicalModelGenerator.makeVolcTopo
2324
GeophysicalModelGenerator.ConstantTemp
2425
GeophysicalModelGenerator.LinearTemp
2526
GeophysicalModelGenerator.HalfspaceCoolingTemp

‎src/LaMEM_io.jl‎

Lines changed: 204 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
using Base: Int64, Float64, NamedTuple
22
using Printf
33
using Glob
4+
using Interpolations
45

56
# LaMEM I/O
67
#
@@ -15,34 +16,42 @@ export GetProcessorPartitioning, ReadData_VTR, ReadData_PVTR, CreatePartitioning
1516
Structure that holds information about the LaMEM grid (usually read from an input file).
1617
"""
1718
struct LaMEM_grid <: AbstractGeneralGrid
19+
# number of markers per element
1820
nmark_x :: Int64
1921
nmark_y :: Int64
2022
nmark_z :: Int64
23+
# total number of markers
2124
nump_x :: Int64
2225
nump_y :: Int64
2326
nump_z :: Int64
24-
27+
# total number of elements in grid
2528
nel_x :: Int64
2629
nel_y :: Int64
2730
nel_z :: Int64
28-
31+
# exent of the grid
2932
W :: Float64
3033
L :: Float64
3134
H :: Float64
32-
35+
# start and end coordinates of grid segments
3336
coord_x
3437
coord_y
3538
coord_z
36-
37-
x1D_c
38-
y1D_c
39-
z1D_c
40-
39+
# 1D vectors with marker coordinates
40+
x_vec
41+
y_vec
42+
z_vec
43+
# grid with marker coordinates
4144
X
4245
Y
4346
Z
44-
45-
47+
# 1D vectors with node coordinates
48+
xn_vec
49+
yn_vec
50+
zn_vec
51+
# grid with node coordinates
52+
Xn
53+
Yn
54+
Zn
4655
end
4756

4857
"""
@@ -141,61 +150,186 @@ z ϵ [-2.0 : 0.0]
141150
"""
142151
function ReadLaMEM_InputFile(file)
143152

144-
nmark_x = ParseValue_LaMEM_InputFile(file,"nmark_x",Int64)
145-
nmark_y = ParseValue_LaMEM_InputFile(file,"nmark_y",Int64)
146-
nmark_z = ParseValue_LaMEM_InputFile(file,"nmark_z",Int64)
153+
# read information from file
154+
nmark_x = ParseValue_LaMEM_InputFile(file,"nmark_x",Int64);
155+
nmark_y = ParseValue_LaMEM_InputFile(file,"nmark_y",Int64);
156+
nmark_z = ParseValue_LaMEM_InputFile(file,"nmark_z",Int64);
157+
158+
nel_x = ParseValue_LaMEM_InputFile(file,"nel_x",Int64);
159+
nel_y = ParseValue_LaMEM_InputFile(file,"nel_y",Int64);
160+
nel_z = ParseValue_LaMEM_InputFile(file,"nel_z",Int64);
161+
162+
coord_x = ParseValue_LaMEM_InputFile(file,"coord_x",Float64);
163+
coord_y = ParseValue_LaMEM_InputFile(file,"coord_y",Float64);
164+
coord_z = ParseValue_LaMEM_InputFile(file,"coord_z",Float64);
165+
166+
# compute infromation from file
167+
W = coord_x[end]-coord_x[1];
168+
L = coord_y[end]-coord_y[1];
169+
H = coord_z[end]-coord_z[1];
170+
171+
nel_x_tot = sum(nel_x);
172+
nel_y_tot = sum(nel_y);
173+
nel_z_tot = sum(nel_z);
174+
175+
nump_x = nel_x_tot*nmark_x;
176+
nump_y = nel_y_tot*nmark_y;
177+
nump_z = nel_z_tot*nmark_z;
178+
179+
# uniform spacing
180+
if (length(coord_x)==2) && (length(coord_y)==2) && (length(coord_z)==2)
181+
# node spacing
182+
dx = W / nel_x_tot;
183+
dy = L / nel_y_tot;
184+
dz = H / nel_z_tot;
185+
186+
# node coordinate vectors
187+
xn = coord_x[1] : dx : coord_x[end];
188+
yn = coord_y[1] : dy : coord_y[end];
189+
zn = coord_z[1] : dz : coord_z[end];
190+
191+
# node grid
192+
Xn,Yn,Zn = XYZGrid(xn, yn, zn);
193+
194+
# marker spacing
195+
dx = W / nump_x;
196+
dy = L / nump_y;
197+
dz = H / nump_z;
198+
199+
# marker coordinate vectors
200+
x = coord_x[1]+dx/2 : dx : coord_x[end]-dx/2;
201+
y = coord_y[1]+dy/2 : dy : coord_y[end]-dy/2;
202+
z = coord_z[1]+dz/2 : dz : coord_z[end]-dz/2;
203+
204+
# marker grid
205+
X,Y,Z = XYZGrid(x, y, z);
206+
207+
# non-uniform spacing
208+
else
209+
# read missing information
210+
nseg_x = ParseValue_LaMEM_InputFile(file,"nseg_x",Int64);
211+
nseg_y = ParseValue_LaMEM_InputFile(file,"nseg_y",Int64);
212+
nseg_z = ParseValue_LaMEM_InputFile(file,"nseg_z",Int64);
213+
214+
bias_x = ParseValue_LaMEM_InputFile(file,"bias_x",Float64);
215+
bias_y = ParseValue_LaMEM_InputFile(file,"bias_y",Float64);
216+
bias_z = ParseValue_LaMEM_InputFile(file,"bias_z",Float64);
217+
218+
## Nodes
219+
# make coordinate vectors
220+
xn = make1DCoords(nseg_x, nel_x, coord_x, bias_x);
221+
yn = make1DCoords(nseg_y, nel_y, coord_y, bias_y);
222+
zn = make1DCoords(nseg_z, nel_z, coord_z, bias_z);
223+
224+
# node grid
225+
Xn,Yn,Zn = XYZGrid(xn, yn, zn);
226+
227+
## Markers
228+
# make marker coordinate vectors
229+
x = make1DMarkerCoords(xn, nmark_x);
230+
y = make1DMarkerCoords(yn, nmark_y);
231+
z = make1DMarkerCoords(zn, nmark_z);
232+
233+
# marker grid
234+
X, Y, Z = XYZGrid(x, y, z);
235+
end
147236

148-
nel_x = ParseValue_LaMEM_InputFile(file,"nel_x",Int64)
149-
nel_y = ParseValue_LaMEM_InputFile(file,"nel_y",Int64)
150-
nel_z = ParseValue_LaMEM_InputFile(file,"nel_z",Int64)
237+
# finish Grid
238+
Grid = LaMEM_grid( nmark_x, nmark_y, nmark_z,
239+
nump_x, nump_y, nump_z,
240+
nel_x_tot, nel_y_tot, nel_z_tot,
241+
W, L, H,
242+
coord_x, coord_y, coord_z,
243+
x, y, z,
244+
X, Y, Z,
245+
xn, yn, zn,
246+
Xn, Yn, Zn);
151247

152-
coord_x = ParseValue_LaMEM_InputFile(file,"coord_x",Float64)
153-
coord_y = ParseValue_LaMEM_InputFile(file,"coord_y",Float64)
154-
coord_z = ParseValue_LaMEM_InputFile(file,"coord_z",Float64)
248+
return Grid
249+
end
155250

156-
if (length(coord_x)>2) || (length(coord_y)>2) || (length(coord_z)>2)
157-
error("Routine currently not working for variable grid spacing")
251+
function make1DMarkerCoords(xn::Array{Float64, 1}, nmark::Int64)
252+
# preallocate
253+
nel = length(xn) - 1
254+
nump = nel * nmark;
255+
x = zeros(Float64, nump);
256+
257+
# compute coordinates
258+
for i = 1 : nel
259+
# start of cell
260+
x0 = xn[i];
261+
# markers spacing inside cell
262+
dx = (xn[i+1] - x0) / nmark;
263+
264+
# compute position
265+
for j = 1 : nmark
266+
x[nmark*i-(nmark-j)] = x0 + dx/2 + (j-1)*dx;
267+
end
158268
end
159269

160-
W = coord_x[end]-coord_x[1];
161-
L = coord_y[end]-coord_y[1];
162-
H = coord_z[end]-coord_z[1];
270+
return x
271+
end
272+
273+
function make1DCoords(nseg::Int64, nel, coord::Array{Float64, 1}, bias)
274+
# preallocate
275+
nel_tot = sum(nel);
276+
x = zeros(Float64, nel_tot+1);
163277

164-
nump_x = nel_x*nmark_x;
165-
nump_y = nel_y*nmark_y;
166-
nump_z = nel_z*nmark_z;
278+
for i = 1 : nseg
279+
# indices of this segment in the coordinate vector
280+
if i == 1
281+
indE = nel[1] + 1
282+
else
283+
indE = sum(nel[1:i]) + 1;
284+
end
285+
indS = indE - nel[i];
167286

168-
dx = W/nump_x;
169-
dy = L/nump_y;
170-
dz = H/nump_z;
287+
# compute coordinates
288+
x[indS:indE] = makeCoordSegment(coord[i], coord[i+1], nel[i], bias[i]);
289+
end
171290

172-
# these lines should be replaced with a separate routine for variable spacing
173-
x = coord_x[1]+dx/2: dx : coord_x[end]-dx/2;
174-
y = coord_y[1]+dy/2: dy : coord_y[end]-dy/2;
175-
z = coord_z[1]+dz/2: dz : coord_z[end]-dz/2;
291+
return x
292+
end
176293

177-
X,Y,Z = XYZGrid(x,y,z); # create 3D grid using regular spacing
178-
179-
Grid = LaMEM_grid( nmark_x, nmark_y, nmark_z,
180-
nump_x, nump_y, nump_z,
181-
nel_x, nel_y, nel_z,
182-
W, L, H,
183-
coord_x, coord_y, coord_z,
184-
x, y, z,
185-
X, Y, Z);
294+
function makeCoordSegment(xStart::Float64, xEnd::Float64, numCells::Int64, bias::Float64)
295+
# average cell size
296+
avgSize = (xEnd - xStart) / numCells;
186297

187-
return Grid
298+
# uniform case
299+
if bias == 1.0
300+
x = Array(xStart : avgSize : xEnd);
301+
# non-uniform case
302+
else
303+
x = zeros(Float64, numCells+1)
304+
# cell size limits
305+
begSize = 2.0 * avgSize / (1.0 + bias);
306+
endSize = bias * begSize;
307+
308+
# cell size increment (negative for bias < 1)
309+
dx = (endSize - begSize) / (numCells - 1);
310+
311+
# generate coordinates
312+
x[1] = xStart;
313+
for i = 2 : numCells + 1
314+
x[i] = x[i-1] + begSize + (i-2)*dx;
315+
end
316+
317+
# overwrite last coordinate
318+
x[end] = xEnd;
319+
end
320+
321+
return x
188322
end
189323

190324
# Print an overview of the LaMEM Grid struct:
191325
function Base.show(io::IO, d::LaMEM_grid)
192326
println(io,"LaMEM Grid: ")
193327
println(io," nel : ($(d.nel_x), $(d.nel_y), $(d.nel_z))")
194328
println(io," marker/cell : ($(d.nmark_x), $(d.nmark_y), $(d.nmark_z))")
195-
println(io," markers : ($(d.nump_x), $(d.nump_x), $(d.nump_x))")
196-
println(io," x ϵ [$(d.coord_x[1]) : $(d.coord_x[2])]")
197-
println(io," y ϵ [$(d.coord_y[1]) : $(d.coord_y[2])]")
198-
println(io," z ϵ [$(d.coord_z[1]) : $(d.coord_z[2])]")
329+
println(io," markers : ($(d.nump_x), $(d.nump_y), $(d.nump_z))")
330+
println(io," x ϵ [$(d.coord_x[1]) : $(d.coord_x[end])]")
331+
println(io," y ϵ [$(d.coord_y[1]) : $(d.coord_y[end])]")
332+
println(io," z ϵ [$(d.coord_z[1]) : $(d.coord_z[end])]")
199333
end
200334

201335
"""
@@ -759,14 +893,30 @@ function Save_LaMEMTopography(Topo::CartData, filename::String)
759893
error("The topography `CartData` structure requires a field :Topography")
760894
end
761895

762-
# Code the topograhic data into a vector
896+
# get grid properties
763897
nx = Float64(size(Topo.fields.Topography,1));
764898
ny = Float64(size(Topo.fields.Topography,2));
765-
x0 = ustrip(Topo.x.val[1,1,1])
766-
y0 = ustrip(Topo.y.val[1,1,1])
767-
dx = ustrip(Topo.x.val[2,2,1]) - x0
768-
dy = ustrip(Topo.y.val[2,2,1]) - y0
769-
Topo_vec = [ nx;ny;x0;y0;dx;dy; ustrip.(Topo.fields.Topography[:])]
899+
x0 = ustrip(Topo.x.val[1,1,1]);
900+
y0 = ustrip(Topo.y.val[1,1,1]);
901+
902+
# LaMEM wants a uniform grid, so interpolate if necessary
903+
if length(unique(trunc.(diff(Topo.x.val[:,1,1]), digits=8))) > 1 || length(unique(trunc.(diff(Topo.y.val[1,:,1]), digits=8))) > 1
904+
x1 = ustrip(Topo.x.val[end,1,1]);
905+
y1 = ustrip(Topo.y.val[1,end,1]);
906+
dx = (x1-x0) / (nx-1);
907+
dy = (y1-y0) / (ny-1);
908+
909+
itp = LinearInterpolation((Topo.x.val[:,1,1], Topo.y.val[1,:,1]), ustrip.(Topo.fields.Topography[:,:,1]));
910+
Topo_itp = [itp(x,y) for x in x0:dx:x1, y in y0:dy:y1];
911+
912+
# Code the topograhic data into a vector
913+
Topo_vec = [ nx;ny;x0;y0;dx;dy; Topo_itp[:]]
914+
else
915+
dx = ustrip(Topo.x.val[2,2,1]) - x0
916+
dy = ustrip(Topo.y.val[2,2,1]) - y0
917+
# Code the topograhic data into a vector
918+
Topo_vec = [ nx;ny;x0;y0;dx;dy; ustrip.(Topo.fields.Topography[:])]
919+
end
770920

771921
# Write as PetscBinary file
772922
PetscBinaryWrite_Vec(filename, Topo_vec)

0 commit comments

Comments
 (0)