From 065c11fd11bc76ada69367134b0d000482cc40d0 Mon Sep 17 00:00:00 2001 From: Seph Date: Fri, 10 Mar 2023 16:50:02 +0100 Subject: [PATCH 01/18] progress on change vtk reader from python to Julia --- src/read_timestep.jl | 253 ++++++++++-------------------------------- test/read_timestep.jl | 8 +- 2 files changed, 62 insertions(+), 199 deletions(-) diff --git a/src/read_timestep.jl b/src/read_timestep.jl index 1c48f438..4f677f5d 100644 --- a/src/read_timestep.jl +++ b/src/read_timestep.jl @@ -4,27 +4,11 @@ # Make these routines easily available outside the module: using GeophysicalModelGenerator: CartData, XYZGrid -using Glob +using Glob, ReadVTK -export Read_VTR_File, field_names, readPVD, Read_VTU_File - - -function vtkXMLPRectilinearGridReader(FileName) - reader1 = pyvtk.vtkXMLPRectilinearGridReader() - reader1.SetFileName(FileName) - reader1.Update() - data = reader1.GetOutput() - return data -end - -function vtkXMLPUnstructuredGridReader(FileName) - reader1 = pyvtk.vtkXMLPUnstructuredGridReader() - reader1.SetFileName(FileName) - reader1.Update() - data = reader1.GetOutput() - return data - end +export Read_LaMEM_PVTR_File, field_names, readPVD, Read_VTU_File + """ output, isCell = ReadField_3D_pVTR(data, FieldName::String) @@ -39,122 +23,25 @@ Output: `data_field` is a tuple of size 1, or 3 depending on whether it is a scalar or vector field """ -function ReadField_3D_pVTR(data, FieldName) - - x = data.GetXCoordinates(); - y = data.GetYCoordinates(); - z = data.GetZCoordinates(); - - nx = PythonCall.pyconvert(Int, x.GetSize()); - ny = PythonCall.pyconvert(Int, y.GetSize()); - nz = PythonCall.pyconvert(Int, z.GetSize()); +function ReadField_3D_pVTR(pvtk, FieldName) isCell = false; - # First assume they are point data - data_f = data.GetPointData().GetArray(FieldName) - if PythonCall.pyisnone(data_f) - data_f = data.GetCellData().GetArray(FieldName) # Try Cell Data - end - numData = PythonCall.pyconvert(Int,data_f.GetDataSize()) - data_Field = PythonCall.pyconvert(Array, data_f); - if size(data_Field,1) != nx*ny*nz - isCell = true; - nx = nx-1; - ny = ny-1; - nz = nz-1; - end + # first try to get point data + data_f = get_point_data(pvtk) + # if empty then load cell data + if isempty(keys(data_f)) + data_f = get_cell_data(pvtk) + data_Field = get_data_reshaped(data_f[FieldName], cell_data=true) - if size(data_Field,2) == 1 - data_Field = reshape(data_Field ,(nx,ny,nz)) if typeof(data_Field[1])==UInt8 data_Field = Int64.(data_Field) end - - data_Tuple = (data_Field, ) - elseif size(data_Field,2) == 3 - Vx = reshape(data_Field[:,1],(nx,ny,nz)); - Vy = reshape(data_Field[:,2],(nx,ny,nz)); - Vz = reshape(data_Field[:,3],(nx,ny,nz)); - data_Tuple = ((Vx,Vy,Vz),) - elseif size(data_Field,2) == 9 - xx = reshape(data_Field[:,1],(nx,ny,nz)); - xy = reshape(data_Field[:,2],(nx,ny,nz)); - xz = reshape(data_Field[:,3],(nx,ny,nz)); - yx = reshape(data_Field[:,4],(nx,ny,nz)); - yy = reshape(data_Field[:,5],(nx,ny,nz)); - yz = reshape(data_Field[:,6],(nx,ny,nz)); - zx = reshape(data_Field[:,7],(nx,ny,nz)); - zy = reshape(data_Field[:,8],(nx,ny,nz)); - zz = reshape(data_Field[:,9],(nx,ny,nz)); - - data_Tuple = ((xx,xy,xz,yx,yy,yz,zx,zy,zz),) + data_Tuple = (data_Field,) + isCell = true; else - error("Not yet implemented for this size $(size(data_Field,2))") - end - name = filter(x -> !isspace(x), FieldName) # remove white spaces - - id = findfirst("[", name) - if !isnothing(id) - name = name[1:id[1]-1] # strip out "[" signs + data_Tuple = (get_data_reshaped(data_f[FieldName]),) end - data_out = NamedTuple{(Symbol(name),)}(data_Tuple,); - - return data_out, isCell -end - -""" - output, isCell = ReadField_3D_pVTU(data, FieldName::String) - -Extracts a 3D data field from a pVTU data structure `data` -Input: -- `data`: Data structure obtained with Read_VTR_File -- `FieldName`: Exact name of the field as specified in the *.vtr file - -Output: -- `data_field`: Array with data, `data_field` is a tuple of size 1, 3 or 9 depending on whether it is a scalar, vector or tensor field - -""" -function ReadField_3D_pVTU(data, FieldName) - - n = PythonCall.pyconvert(Int, data.GetNumberOfPoints()); - coords = PythonCall.pyconvert(Array,data.GetPoints().GetData()); # coordinates of points - - # First assume they are point data - data_f = data.GetPointData().GetArray(FieldName) - if PythonCall.pyisnone(data_f) - data_f = data.GetCellData().GetArray(FieldName) # Try Cell Data - end - numData = PythonCall.pyconvert(Int,data_f.GetDataSize()) - data_Field = PythonCall.pyconvert(Array, data_f); - - if size(data_Field,2) == 1 - data_Field = data_Field; - if typeof(data_Field[1])==UInt8 - data_Field = Int64.(data_Field) - end - - data_Tuple = (data_Field, ) - elseif size(data_Field,2) == 3 - Vx = data_Field[:,1] - Vy = data_Field[:,2] - Vz = data_Field[:,3] - data_Tuple = ((Vx,Vy,Vz),) - elseif size(data_Field,2) == 9 - xx = data_Field[:,1] - xy = data_Field[:,2] - xz = data_Field[:,3] - yx = data_Field[:,4] - yy = data_Field[:,5] - yz = data_Field[:,6] - zx = data_Field[:,7] - zy = data_Field[:,8] - zz = data_Field[:,9] - - data_Tuple = ((xx,xy,xz,yx,yy,yz,zx,zy,zz),) - else - error("Not yet implemented for this size $(size(data_Field,2))") - end name = filter(x -> !isspace(x), FieldName) # remove white spaces id = findfirst("[", name) @@ -163,60 +50,12 @@ function ReadField_3D_pVTU(data, FieldName) end data_out = NamedTuple{(Symbol(name),)}(data_Tuple,); - return data_out -end - - - -""" - names = field_names(data) - -Returns the names of all fields stored in the vtr data structure `data`. -""" -function field_names(data) - names = []; - - # Get Names of PointData arrays - pdata = data.GetPointData() - num = PythonCall.pyconvert(Int,pdata.GetNumberOfArrays()) - for i=1:num - names = [names; PythonCall.pyconvert(String,pdata.GetArrayName(i-1))] - end - - # Get Names of CellData arrays - cdata = data.GetCellData() - num = PythonCall.pyconvert(Int,cdata.GetNumberOfArrays()) - for i=1:num - names = [names; PythonCall.pyconvert(String,cdata.GetArrayName(i-1))] - end - - return names; -end - -""" - names = field_names(DirName, FileName) - -Returns the names of all fields stored in the vtr file in the directory `Directory` and file `FileName`. -""" -function field_names(DirName, FileName) - CurDir = pwd(); - - # change to directory - cd(DirName) - - # read data from parallel rectilinear grid - data = vtkXMLPRectilinearGridReader(FileName); # this is how it should be done within modules - cd(CurDir) - - # Fields stored in this data file: - names = field_names(data); - - return names; + return data_out, isCell end """ - data_output = Read_VTR_File(DirName, FileName; field=nothing) + data_output = Read_LaMEM_PVTR_File(DirName, FileName; field=nothing) Reads a 3D LaMEM timestep from VTR file `FileName`, located in directory `DirName`. By default, it will read all fields. If you want you can only read a specific `field`. See the function `fieldnames` to get a list with all available fields in the file. @@ -224,24 +63,30 @@ By default, it will read all fields. If you want you can only read a specific `f It will return `data_output` which is a `CartData` output structure. """ -function Read_VTR_File(DirName, FileName; field=nothing) +function Read_LaMEM_PVTR_File(DirName, FileName; field=nothing) CurDir = pwd(); # change to directory cd(DirName) # read data from parallel rectilinear grid - data = vtkXMLPRectilinearGridReader(FileName); # this is how it should be done within modules + pvtk = PVTKFile(FileName) + cd(CurDir) - # Fields stored in this data file: - names = field_names(data); + # Fields stored in this data file: + names = keys(get_point_data(pvtk)) + if isempty(names) + names = keys(get_cell_data(pvtk)) + end + isCell = true if isnothing(field) # read all data in the file data_fields = NamedTuple(); + for FieldName in names - dat, isCell = ReadField_3D_pVTR(data, FieldName); + dat, isCell = ReadField_3D_pVTR(pvtk, FieldName); data_fields = merge(data_fields,dat) end @@ -256,10 +101,13 @@ function Read_VTR_File(DirName, FileName; field=nothing) end end + coords_read = get_coordinates(pvtk) + # Read coordinates - x = PythonCall.pyconvert(Array,data.GetXCoordinates()); - y = PythonCall.pyconvert(Array,data.GetYCoordinates()); - z = PythonCall.pyconvert(Array,data.GetZCoordinates()); + x = coords_read[1] + y = coords_read[2] + z = coords_read[3] + if isCell # In case we have cell data , coordinates are center of cell x = (x[1:end-1] + x[2:end])/2 @@ -272,6 +120,7 @@ function Read_VTR_File(DirName, FileName; field=nothing) return data_output end + """ FileNames, Time = readPVD(FileName::String) @@ -315,17 +164,22 @@ function Read_VTU_File(DirName, FileName; field=nothing) cd(DirName) # read data from parallel rectilinear grid - data = vtkXMLPUnstructuredGridReader(FileName); # this is how it should be done within modules + pvtu = PVTKFile(FileName) + cd(CurDir) - # Fields stored in this data file: - names = field_names(data); + # Fields stored in this data file: + names = keys(get_point_data(pvtu)) + if isempty(names) + names = keys(get_cell_data(pvtu)) + end + isCell = true if isnothing(field) # read all data in the file data_fields = NamedTuple(); for FieldName in names - dat = ReadField_3D_pVTU(data, FieldName); + dat = ReadField_3D_pVTU(pvtu, FieldName); data_fields = merge(data_fields,dat) end @@ -340,12 +194,21 @@ function Read_VTU_File(DirName, FileName; field=nothing) end end + coords_read = get_coordinates(pvtu) + # Read coordinates - coords = PythonCall.pyconvert(Array,data.GetPoints().GetData()); # coordinates of points - x = coords[:,1]; - y = coords[:,2]; - z = coords[:,3]; - - data_output = CartData(x,y,z, data_fields) - return data_output + x = coords_read[1] + y = coords_read[2] + z = coords_read[3] + + if isCell + # In case we have cell data , coordinates are center of cell + x = (x[1:end-1] + x[2:end])/2 + y = (y[1:end-1] + y[2:end])/2 + z = (z[1:end-1] + z[2:end])/2 + end + + X,Y,Z = XYZGrid(x,y,z) + data_output = CartData(X,Y,Z, data_fields) + return data_output end \ No newline at end of file diff --git a/test/read_timestep.jl b/test/read_timestep.jl index 7fc5f1be..58819aae 100644 --- a/test/read_timestep.jl +++ b/test/read_timestep.jl @@ -7,20 +7,20 @@ using LaMEM FileName="FB_multigrid.pvtr" DirName = "Timestep_00000001_6.72970343e+00" - data = Read_VTR_File(DirName, FileName) + data = Read_LaMEM_PVTR_File(DirName, FileName) @test sum(data.fields.phase) ≈ 736.36414f0 - @test sum(data.fields.strain_rate[1]) ≈ -0.019376338f0 + @test sum(data.fields.strain_rate[1,:,:,:]) ≈ -0.019376338f0 fields = field_names(DirName, FileName) # with cell-data FileName="FB_multigrid_phase.pvtr" DirName = "Timestep_00000001_6.72970343e+00" - data = Read_VTR_File(DirName, FileName) + data = Read_LaMEM_PVTR_File(DirName, FileName) @test sum(data.fields.phase) == 19822 # read subduction setup - data = Read_VTR_File("Timestep_00000001_5.50000000e-02", "Subduction2D_FreeSlip_direct.pvtr") + data = Read_LaMEM_PVTR_File("Timestep_00000001_5.50000000e-02", "Subduction2D_FreeSlip_direct.pvtr") @test sum(data.fields.density) ≈ 1.60555f8 # Read PVD file From 4b0d1b0d14857eb93ea3ed892c1326c58f863298 Mon Sep 17 00:00:00 2001 From: Seph Date: Fri, 10 Mar 2023 21:42:00 +0100 Subject: [PATCH 02/18] adding pvts reader --- src/read_timestep.jl | 148 ++++++++++++- .../Subduction2D_FreeSurface_DirectSolver.dat | 201 ++++++++++++++++++ test/read_timestep.jl | 12 +- test/runLaMEM.jl | 6 + 4 files changed, 357 insertions(+), 10 deletions(-) create mode 100644 test/input_files/Subduction2D_FreeSurface_DirectSolver.dat diff --git a/src/read_timestep.jl b/src/read_timestep.jl index 4f677f5d..27b61770 100644 --- a/src/read_timestep.jl +++ b/src/read_timestep.jl @@ -6,7 +6,7 @@ using GeophysicalModelGenerator: CartData, XYZGrid using Glob, ReadVTK -export Read_LaMEM_PVTR_File, field_names, readPVD, Read_VTU_File +export Read_LaMEM_PVTR_File, Read_LaMEM_PVTS_File, field_names, readPVD, Read_LaMEM_PVTU_File """ @@ -54,6 +54,80 @@ function ReadField_3D_pVTR(pvtk, FieldName) end + +function ReadField_3D_pVTS(pvts, FieldName) + isCell = false; + + # first try to get point data + data_f = get_point_data(pvts) + # if empty then load cell data + if isempty(keys(data_f)) + data_f = get_cell_data(pvts) + data_Field = get_data_reshaped(data_f[FieldName], cell_data=true)[:,:,:,1] + + if typeof(data_Field[1])==UInt8 + data_Field = Int64.(data_Field) + end + data_Tuple = (data_Field,) + isCell = true; + else + data_Tuple = (get_data_reshaped(data_f[FieldName])[:,:,:,1],) + end + + name = filter(x -> !isspace(x), FieldName) # remove white spaces + + id = findfirst("[", name) + if !isnothing(id) + name = name[1:id[1]-1] # strip out "[" signs + end + data_out = NamedTuple{(Symbol(name),)}(data_Tuple,); + + return data_out, isCell +end + + +""" + output, isCell = ReadField_3D_pVTU(data, FieldName::String) +Extracts a 3D data field from a pVTU data structure `data` +Input: +- `data`: Data structure obtained with Read_VTR_File +- `FieldName`: Exact name of the field as specified in the *.vtr file + +Output: +- `data_field`: Array with data, `data_field` is a tuple of size 1, 3 or 9 depending on whether it is a scalar, vector or tensor field + +""" +function ReadField_3D_pVTU(pvtu, FieldName) + + # first try to get point data + data_f = get_point_data(pvtu) + # if empty then load cell data + if isempty(keys(data_f)) + data_f = get_cell_data(pvtk) + data_Field = get_data(data_f[FieldName], cell_data=true)[1] + + if typeof(data_Field[1])==UInt8 + data_Field = Int64.(data_Field) + end + data_Tuple = (data_Field,) + + else + data_Tuple = (get_data(data_f[FieldName])[1],) + end + + name = filter(x -> !isspace(x), FieldName) # remove white spaces + + id = findfirst("[", name) + if !isnothing(id) + name = name[1:id[1]-1] # strip out "[" signs + end + data_out = NamedTuple{(Symbol(name),)}(data_Tuple,); + + return data_out +end + + + """ data_output = Read_LaMEM_PVTR_File(DirName, FileName; field=nothing) @@ -149,7 +223,7 @@ end """ - data_output = Read_VTU_File(DirName, FileName; field=nothing) + data_output = Read_LaMEM_PVTU_File(DirName, FileName; field=nothing) Reads a 3D LaMEM timestep from VTU file `FileName`, located in directory `DirName`. Typically this is done to read passive tracers back into julia. By default, it will read all fields. If you want you can only read a specific `field`. See the function `fieldnames` to get a list with all available fields in the file. @@ -157,7 +231,7 @@ By default, it will read all fields. If you want you can only read a specific `f It will return `data_output` which is a `CartData` output structure. """ -function Read_VTU_File(DirName, FileName; field=nothing) +function Read_LaMEM_PVTU_File(DirName, FileName; field=nothing) CurDir = pwd(); # change to directory @@ -194,12 +268,70 @@ function Read_VTU_File(DirName, FileName; field=nothing) end end - coords_read = get_coordinates(pvtu) + points = get_points(pvtu) # Read coordinates - x = coords_read[1] - y = coords_read[2] - z = coords_read[3] + x = points[1][1,:] + y = points[1][2,:] + z = points[1][3,:] + + data_output = CartData(x,y,z, data_fields) + return data_output +end + + +""" + data_output = Read_LaMEM_PVTU_File(DirName, FileName; field=nothing) + +Reads a 3D LaMEM timestep from VTU file `FileName`, located in directory `DirName`. Typically this is done to read passive tracers back into julia. +By default, it will read all fields. If you want you can only read a specific `field`. See the function `fieldnames` to get a list with all available fields in the file. + +It will return `data_output` which is a `CartData` output structure. + +""" +function Read_LaMEM_PVTS_File(DirName, FileName; field=nothing) + CurDir = pwd(); + + # change to directory + cd(DirName) + + # read data from parallel rectilinear grid + pvts = PVTKFile(FileName) + + cd(CurDir) + + # Fields stored in this data file: + name = keys(get_point_data(pvts)) + if isempty(name) + name = keys(get_cell_data(pvts)) + end + + isCell = true + if isnothing(field) + # read all data in the file + data_fields = NamedTuple(); + for FieldName in name + dat, isCell = ReadField_3D_pVTS(pvts, FieldName); + data_fields = merge(data_fields,dat) + end + + else + # read just a single data set + ind = findall(contains.(field, name)) + # check that it exists + if isempty(ind) + error("the field $field does not exist in the data file") + else + data_fields, isCell = ReadField_3D_pVTS(data, field) + end + end + + coords_read = get_coordinates(pvts) + + # Read coordinates + x = coords_read[1][:,1] + y = coords_read[2][:,1] + z = coords_read[3][:,1] if isCell # In case we have cell data , coordinates are center of cell @@ -210,5 +342,5 @@ function Read_VTU_File(DirName, FileName; field=nothing) X,Y,Z = XYZGrid(x,y,z) data_output = CartData(X,Y,Z, data_fields) - return data_output + return data_output end \ No newline at end of file diff --git a/test/input_files/Subduction2D_FreeSurface_DirectSolver.dat b/test/input_files/Subduction2D_FreeSurface_DirectSolver.dat new file mode 100644 index 00000000..6b1d1d07 --- /dev/null +++ b/test/input_files/Subduction2D_FreeSurface_DirectSolver.dat @@ -0,0 +1,201 @@ +# This is a simple 2D subduction setup with GEO units and a free surface upper boundary +# The model is linear viscous and the geometry is created with the build-in geometry objects +# The full simulation runs for 150 Myrs, and takes about 500 timesteps + + +#=============================================================================== +# Scaling +#=============================================================================== + + units = geo # geological units + + unit_temperature = 1000 + unit_length = 1e3 + unit_viscosity = 1e20 + unit_stress = 1e9 + +#=============================================================================== +# Time stepping parameters +#=============================================================================== + + time_end = 100 # simulation end time + dt = 0.1 # initial time step + dt_min = 1e-5 # minimum time step (declare divergence if lower value is attempted) + dt_max = 100 # maximum time step + inc_dt = 0.1 # time step increment per time step (fraction of unit) + CFL = 0.5 # CFL (Courant-Friedrichs-Lewy) criterion + CFLMAX = 0.5 # CFL criterion for elasticity + nstep_max = 1000 # maximum allowed number of steps (lower bound: time_end/dt_max) + nstep_out = 5 # save output every n steps + nstep_rdb = 100 # save restart database every n steps + + +#=============================================================================== +# Grid & discretization parameters +#=============================================================================== + +# Number of cells for all segments + nel_x = 256 + nel_y = 2 + nel_z = 64 + +# Coordinates of all segments (including start and end points) + + coord_x = -1500 1500 + coord_y = -10 10 + coord_z = -660 40 + +#=============================================================================== +# Free surface +#=============================================================================== + surf_use = 1 # free surface activation flag + surf_corr_phase = 1 # air phase ratio correction flag (due to surface position) + surf_level = 0 # initial level + surf_air_phase = 2 # phase ID of sticky air layer + surf_max_angle = 10.0 # maximum angle with horizon (smoothed if larger) + +#=============================================================================== +# Boundary conditions +#=============================================================================== +# No-slip boundary flag mask (left right front back bottom top) + + noslip = 0 0 0 0 1 0 + + +#=============================================================================== +# Solution parameters & controls +#=============================================================================== + + gravity = 0.0 0.0 -9.81 # gravity vector + FSSA = 1.0 # free surface stabilization parameter [0 - 1] + init_guess = 0 # initial guess flag + eta_min = 1e18 # viscosity upper bound + eta_ref = 1e20 # reference viscosity for initial guess + eta_max = 1e25 # viscosity lower limit + + +#=============================================================================== +# Solver options +#=============================================================================== + SolverType = direct # solver [direct or multigrid] + DirectSolver = mumps # mumps/superlu_dist/pastix + DirectPenalty = 1e5 + + +#=============================================================================== +# Model setup & advection +#=============================================================================== + + msetup = geom # setup type + nmark_x = 3 # markers per cell in x-direction + nmark_y = 3 # ... y-direction + nmark_z = 3 # ... z-direction + bg_phase = 0 # background phase ID + rand_noise = 1 # random noise flag + advect = rk2 # advection scheme + interp = minmod # velocity interpolation scheme + mark_ctrl = basic # marker control type + nmark_lim = 8 100 # min/max number per cell + + +# Geometric primitives: + # sticky air + + phase = 2 + bounds = -1500 1500 -25 25 0 100 # (left, right, front, back, bottom, top) + + + # horizontal part of slab + + phase = 1 + bounds = 0 1200 -25 25 -80 0 # (left, right, front, back, bottom, top) + + + # Inclined part of slab + + phase = 1 + coord = -250 -25 -230 0 -25 -80 0 25 -80 -250 25 -230 -250 -25 -150 0 -25 0 0 25 0 -250 25 -150 + + +#=============================================================================== +# Output +#=============================================================================== + +# Grid output options (output is always active) + + out_file_name = Subduction2D_FreeSurface_direct # output file name + out_pvd = 1 # activate writing .pvd file + out_phase = 1 + out_density = 1 + out_visc_total = 1 + out_visc_creep = 1 + out_velocity = 1 + out_pressure = 1 + out_eff_press = 1 + out_temperature = 1 + out_dev_stress = 1 + out_j2_dev_stress = 1 + out_strain_rate = 1 + out_j2_strain_rate = 1 + out_yield = 1 + out_plast_strain = 1 + out_plast_dissip = 1 + out_tot_displ = 1 + out_moment_res = 1 + out_cont_res = 1 + +# AVD phase viewer output options (requires activation) + + out_avd = 1 # activate AVD phase output + out_avd_pvd = 1 # activate writing .pvd file + out_avd_ref = 3 # AVD grid refinement factor + +# Free surface output options (can be activated only if surface tracking is enabled) + + out_surf = 1 # activate surface output + out_surf_pvd = 1 # activate writing .pvd file + out_surf_velocity = 1 + out_surf_topography = 1 + out_surf_amplitude = 1 + +#=============================================================================== +# Material phase parameters +#=============================================================================== + + # Define properties of mantle + + ID = 0 # phase id + rho = 3200 # density + eta = 1e21 # viscosity + + + # Define properties of slab + + ID = 1 # phase id + rho = 3300 # density + eta = 1e23 # viscosity + + + # Define properties of sticky air + + ID = 2 # phase id + rho = 1 # density + eta = 1e18 # viscosity + + +#=============================================================================== +# PETSc options +#=============================================================================== + + + + # LINEAR & NONLINEAR SOLVER OPTIONS + -snes_ksp_ew + + + -js_ksp_monitor # display how the inner iterations converge + + + + +#=============================================================================== diff --git a/test/read_timestep.jl b/test/read_timestep.jl index 58819aae..422f61ae 100644 --- a/test/read_timestep.jl +++ b/test/read_timestep.jl @@ -28,8 +28,16 @@ using LaMEM @test Time[2] ≈ 0.055 # Read passive tracers - data = Read_VTU_File("Timestep_00000010_1.09635548e+00", "PlumeLithosphereInteraction_passive_tracers.pvtu") - @test data.z[100] ≈ -298.5178f0 + data = Read_LaMEM_PVTU_File("Timestep_00000010_1.09635548e+00", "PlumeLithosphereInteraction_passive_tracers.pvtu") + @test data.z[100] ≈ -298.4531f0 @test data.fields.Temperature[100] ≈ 1350.0f0 + # Read surface data + data = Read_LaMEM_PVTS_File("Timestep_00000005_7.59607376e-02", "Subduction2D_FreeSurface_direct_surf.pvts") + # @test data.z[100] ≈ -298.4531f0 + # @test data.fields.Temperature[100] ≈ 1350.0f0 + + + end + diff --git a/test/runLaMEM.jl b/test/runLaMEM.jl index 7d1f565e..fb7e8e00 100644 --- a/test/runLaMEM.jl +++ b/test/runLaMEM.jl @@ -29,6 +29,12 @@ using LaMEM out = run_lamem(ParamFile, 1, "-nstep_max 2") # 1 core @test isnothing(out) + # Try free surface + ParamFile="input_files/Subduction2D_FreeSurface_DirectSolver.dat"; + out = run_lamem(ParamFile, 4, "-nstep_max 5") # 4 core + @test isnothing(out) + + if !Sys.iswindows() out = run_lamem(ParamFile, 2, "-nstep_max 2") # 2 cores (mumps) @test isnothing(out) From be7d35bbd45a553c777d7b344e54e5735cbd5c59 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Fri, 10 Mar 2023 22:35:17 +0100 Subject: [PATCH 03/18] remove python; fix vts --- Project.toml | 5 +---- src/LaMEM.jl | 9 +++++---- src/read_timestep.jl | 27 +++++++++++++++------------ 3 files changed, 21 insertions(+), 20 deletions(-) diff --git a/Project.toml b/Project.toml index abd8f84f..8aec278b 100644 --- a/Project.toml +++ b/Project.toml @@ -7,15 +7,12 @@ version = "0.1.6" GeophysicalModelGenerator = "3700c31b-fa53-48a6-808a-ef22d5a84742" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" LaMEM_jll = "15d6fa20-f789-5486-b71b-22b4ac8eb1c1" -MicroMamba = "0b3b1443-0f03-428d-bdfb-f27f9c1191ea" -PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" +ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3" [compat] GeophysicalModelGenerator = "0.4" Glob = "1" LaMEM_jll = "1.2" -MicroMamba = "0.1.9 - 0.1.9" # MicroMamba 0.1.10 appears to be broken -PythonCall = "0.9" julia = "1.6" [extras] diff --git a/src/LaMEM.jl b/src/LaMEM.jl index a579e1ee..f4a07540 100644 --- a/src/LaMEM.jl +++ b/src/LaMEM.jl @@ -1,9 +1,9 @@ module LaMEM using LaMEM_jll -using PythonCall +#using PythonCall using Glob -ENV["JULIA_CONDAPKG_BACKEND"] = "MicroMamba" +#ENV["JULIA_CONDAPKG_BACKEND"] = "MicroMamba" # load the correct MPI const mpiexec = if isdefined(LaMEM_jll,:MPICH_jll) @@ -14,6 +14,7 @@ else nothing end +#= # Reading files back into julia const pyvtk = PythonCall.pynew() @@ -26,7 +27,7 @@ function __init__() # link vtk. Note that all python dependencies are listed in PythonCallDeps.toml PythonCall.pycopy!(pyvtk, pyimport("vtk")) # used to read VTK files end - +=# include("run_lamem.jl") include("run_lamem_save_grid.jl") @@ -35,6 +36,6 @@ include("utils.jl") export run_lamem export run_lamem_save_grid -export pyvtk +#export pyvtk end # module diff --git a/src/read_timestep.jl b/src/read_timestep.jl index 27b61770..6de656d8 100644 --- a/src/read_timestep.jl +++ b/src/read_timestep.jl @@ -8,7 +8,6 @@ using Glob, ReadVTK export Read_LaMEM_PVTR_File, Read_LaMEM_PVTS_File, field_names, readPVD, Read_LaMEM_PVTU_File - """ output, isCell = ReadField_3D_pVTR(data, FieldName::String) @@ -281,9 +280,9 @@ end """ - data_output = Read_LaMEM_PVTU_File(DirName, FileName; field=nothing) + data_output = Read_LaMEM_PVTS_File(DirName, FileName; field=nothing) -Reads a 3D LaMEM timestep from VTU file `FileName`, located in directory `DirName`. Typically this is done to read passive tracers back into julia. +Reads a 3D LaMEM timestep from VTS file `FileName`, located in directory `DirName`. Typically this is done to read passive tracers back into julia. By default, it will read all fields. If you want you can only read a specific `field`. See the function `fieldnames` to get a list with all available fields in the file. It will return `data_output` which is a `CartData` output structure. @@ -306,13 +305,23 @@ function Read_LaMEM_PVTS_File(DirName, FileName; field=nothing) name = keys(get_cell_data(pvts)) end + coords_read = get_coordinates(pvts) + + # Read coordinates + X = coords_read[1] + Y = coords_read[2] + Z = coords_read[3] + isCell = true if isnothing(field) # read all data in the file data_fields = NamedTuple(); for FieldName in name - dat, isCell = ReadField_3D_pVTS(pvts, FieldName); - data_fields = merge(data_fields,dat) + dat, isCell = LaMEM.ReadField_3D_pVTS(pvts, FieldName); + if length(dat[1])>length(X) + dat = NamedTuple{keys(dat)}((ntuple(i->dat[1][i,:,:,:], size(dat[1],1)),)) + end + data_fields = merge(data_fields,dat) end else @@ -326,12 +335,7 @@ function Read_LaMEM_PVTS_File(DirName, FileName; field=nothing) end end - coords_read = get_coordinates(pvts) - - # Read coordinates - x = coords_read[1][:,1] - y = coords_read[2][:,1] - z = coords_read[3][:,1] + if isCell # In case we have cell data , coordinates are center of cell @@ -340,7 +344,6 @@ function Read_LaMEM_PVTS_File(DirName, FileName; field=nothing) z = (z[1:end-1] + z[2:end])/2 end - X,Y,Z = XYZGrid(x,y,z) data_output = CartData(X,Y,Z, data_fields) return data_output end \ No newline at end of file From 3a8f70df8132c901ffd193a3793619fb77f6591a Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Fri, 10 Mar 2023 22:59:01 +0100 Subject: [PATCH 04/18] simplify routine --- src/read_timestep.jl | 58 +++++++++++++------------------------------- 1 file changed, 17 insertions(+), 41 deletions(-) diff --git a/src/read_timestep.jl b/src/read_timestep.jl index 6de656d8..b310aefd 100644 --- a/src/read_timestep.jl +++ b/src/read_timestep.jl @@ -291,57 +291,33 @@ It will return `data_output` which is a `CartData` output structure. function Read_LaMEM_PVTS_File(DirName, FileName; field=nothing) CurDir = pwd(); - # change to directory - cd(DirName) - # read data from parallel rectilinear grid + cd(DirName) pvts = PVTKFile(FileName) - cd(CurDir) # Fields stored in this data file: name = keys(get_point_data(pvts)) - if isempty(name) - name = keys(get_cell_data(pvts)) - end - - coords_read = get_coordinates(pvts) - - # Read coordinates - X = coords_read[1] - Y = coords_read[2] - Z = coords_read[3] - - isCell = true - if isnothing(field) - # read all data in the file - data_fields = NamedTuple(); - for FieldName in name - dat, isCell = LaMEM.ReadField_3D_pVTS(pvts, FieldName); - if length(dat[1])>length(X) - dat = NamedTuple{keys(dat)}((ntuple(i->dat[1][i,:,:,:], size(dat[1],1)),)) - end - data_fields = merge(data_fields,dat) - end - - else - # read just a single data set - ind = findall(contains.(field, name)) - # check that it exists - if isempty(ind) - error("the field $field does not exist in the data file") + if !isnothing(field) + if any(occursin.(name,field)) + name = tuple(field); else - data_fields, isCell = ReadField_3D_pVTS(data, field) + error("the field $field does not exist in the data file which has fields: $(name)") end end - - - if isCell - # In case we have cell data , coordinates are center of cell - x = (x[1:end-1] + x[2:end])/2 - y = (y[1:end-1] + y[2:end])/2 - z = (z[1:end-1] + z[2:end])/2 + # Read coordinates + X,Y,Z = get_coordinates(pvts) + + # read all data in the file + data_fields = NamedTuple(); + for FieldName in name + dat, _ = LaMEM.ReadField_3D_pVTS(pvts, FieldName); + if length(dat[1])>length(X) + dat_tup = ntuple(i->dat[1][i,:,:,:], size(dat[1],1)) + dat = NamedTuple{keys(dat)}(tuple(dat_tup)) + end + data_fields = merge(data_fields,dat) end data_output = CartData(X,Y,Z, data_fields) From 93bb34da0eb35355d24efc92f68e3aff45b6d151 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Fri, 10 Mar 2023 23:05:32 +0100 Subject: [PATCH 05/18] remove python --- src/LaMEM.jl | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/src/LaMEM.jl b/src/LaMEM.jl index f4a07540..c53b5cc7 100644 --- a/src/LaMEM.jl +++ b/src/LaMEM.jl @@ -1,10 +1,7 @@ module LaMEM using LaMEM_jll -#using PythonCall using Glob -#ENV["JULIA_CONDAPKG_BACKEND"] = "MicroMamba" - # load the correct MPI const mpiexec = if isdefined(LaMEM_jll,:MPICH_jll) LaMEM_jll.MPICH_jll.mpiexec() @@ -14,21 +11,6 @@ else nothing end -#= -# Reading files back into julia -const pyvtk = PythonCall.pynew() - -function __init__() - - println("Adding PythonCall dependencies to read LaMEM timesteps") - pth = (@__DIR__)*"/python" # Path where the python routines are - pyimport("sys").path.append(pth) # append path - - # link vtk. Note that all python dependencies are listed in PythonCallDeps.toml - PythonCall.pycopy!(pyvtk, pyimport("vtk")) # used to read VTK files -end -=# - include("run_lamem.jl") include("run_lamem_save_grid.jl") include("read_timestep.jl") @@ -36,6 +18,5 @@ include("utils.jl") export run_lamem export run_lamem_save_grid -#export pyvtk end # module From 080d0fd48f2f9ac89f71f5203f2d9799e0f24b9a Mon Sep 17 00:00:00 2001 From: Seph Date: Sat, 11 Mar 2023 11:21:03 +0100 Subject: [PATCH 06/18] added test for surf file --- test/read_timestep.jl | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/test/read_timestep.jl b/test/read_timestep.jl index 422f61ae..4e3f4de7 100644 --- a/test/read_timestep.jl +++ b/test/read_timestep.jl @@ -34,10 +34,7 @@ using LaMEM # Read surface data data = Read_LaMEM_PVTS_File("Timestep_00000005_7.59607376e-02", "Subduction2D_FreeSurface_direct_surf.pvts") - # @test data.z[100] ≈ -298.4531f0 - # @test data.fields.Temperature[100] ≈ 1350.0f0 - - - + @test data.z[100] ≈ 0.6830405f0 + @test sum(data.fields.topography[:,1,1]) ≈ 1.2634416f0 end From bbb3445d1d1f22bcae21d5a67b0bd8d60214c7bf Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Tue, 14 Mar 2023 10:20:19 +0100 Subject: [PATCH 07/18] update --- src/utils.jl | 1 - test/runtests.jl | 1 - 2 files changed, 2 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 1af9a5a6..62b3c7f7 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -46,7 +46,6 @@ The downloaded `LaMEM` binaries can also be called from outside julia (directly In that case, you will need to set load correct dynamic libraries (such as PETSc) and call the correct binaries. This function shows this for your system. - """ function show_paths_LaMEM() diff --git a/test/runtests.jl b/test/runtests.jl index b0683930..f2f95a2a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,5 @@ using Test using LaMEM -using PythonCall include("runLaMEM.jl") include("read_timestep.jl") From 9a2e8b538c74440b63f857d3a772b2e39f8033f7 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Sat, 8 Apr 2023 19:34:41 +0200 Subject: [PATCH 08/18] changes to put it in modules and make it easier to read a LaMEM timestep --- src/IO.jl | 11 +++++ src/LaMEM.jl | 30 +++++++++--- src/Run.jl | 22 +++++++++ src/read_timestep.jl | 72 +++++++++++++++++++++++++---- src/run_lamem_save_grid.jl | 42 +++++++++-------- src/utils_IO.jl | 80 ++++++++++++++++++++++++++++++++ src/{utils.jl => utils_Run.jl} | 84 +--------------------------------- 7 files changed, 221 insertions(+), 120 deletions(-) create mode 100644 src/IO.jl create mode 100644 src/Run.jl create mode 100644 src/utils_IO.jl rename src/{utils.jl => utils_Run.jl} (53%) diff --git a/src/IO.jl b/src/IO.jl new file mode 100644 index 00000000..e9cd9fdc --- /dev/null +++ b/src/IO.jl @@ -0,0 +1,11 @@ +module IO +# this contains I/O routines of LaMEM, which don't require LaMEM_jll + +include("read_timestep.jl") +export Read_LaMEM_PVTR_File, Read_LaMEM_PVTS_File, field_names, readPVD, Read_LaMEM_PVTU_File + +include("utils_IO.jl") +export clean_directory, changefolder + + +end \ No newline at end of file diff --git a/src/LaMEM.jl b/src/LaMEM.jl index c53b5cc7..328bd47d 100644 --- a/src/LaMEM.jl +++ b/src/LaMEM.jl @@ -1,6 +1,21 @@ module LaMEM + +# Functions to read LaMEM output +include("IO.jl") +using .IO +export Read_LaMEM_PVTR_File, Read_LaMEM_PVTS_File, field_names, readPVD, Read_LaMEM_PVTU_File +export Read_LaMEM_simulation, Read_LaMEM_timestep + +# Functions +include("Run.jl") +using .Run +export run_lamem, run_lamem_save_grid, mpiexec +export remove_popup_messages_mac, show_paths_LaMEM + + + +#= using LaMEM_jll -using Glob # load the correct MPI const mpiexec = if isdefined(LaMEM_jll,:MPICH_jll) @@ -11,12 +26,13 @@ else nothing end -include("run_lamem.jl") -include("run_lamem_save_grid.jl") -include("read_timestep.jl") -include("utils.jl") -export run_lamem -export run_lamem_save_grid + +=# + + +#include("read_timestep.jl") +#include("utils.jl") + end # module diff --git a/src/Run.jl b/src/Run.jl new file mode 100644 index 00000000..81b61775 --- /dev/null +++ b/src/Run.jl @@ -0,0 +1,22 @@ +module Run +# module to run LaMEM_jll +using LaMEM_jll,Glob + +export run_lamem, run_lamem_save_grid, mpiexec +export remove_popup_messages_mac, show_paths_LaMEM + +include("run_lamem.jl") +include("run_lamem_save_grid.jl") +include("utils_Run.jl") + +# load the correct MPI +const mpiexec = if isdefined(LaMEM_jll,:MPICH_jll) + LaMEM_jll.MPICH_jll.mpiexec() +elseif isdefined(LaMEM_jll,:MicrosoftMPI_jll) + LaMEM_jll.MicrosoftMPI_jll.mpiexec() +else + nothing +end + + +end \ No newline at end of file diff --git a/src/read_timestep.jl b/src/read_timestep.jl index b310aefd..bb8fab66 100644 --- a/src/read_timestep.jl +++ b/src/read_timestep.jl @@ -6,7 +6,8 @@ using GeophysicalModelGenerator: CartData, XYZGrid using Glob, ReadVTK -export Read_LaMEM_PVTR_File, Read_LaMEM_PVTS_File, field_names, readPVD, Read_LaMEM_PVTU_File +export Read_LaMEM_PVTR_File, Read_LaMEM_PVTS_File, field_names, Read_LaMEM_PVTU_File +export Read_LaMEM_simulation, Read_LaMEM_timestep """ output, isCell = ReadField_3D_pVTR(data, FieldName::String) @@ -53,7 +54,6 @@ function ReadField_3D_pVTR(pvtk, FieldName) end - function ReadField_3D_pVTS(pvts, FieldName) isCell = false; @@ -136,7 +136,7 @@ By default, it will read all fields. If you want you can only read a specific `f It will return `data_output` which is a `CartData` output structure. """ -function Read_LaMEM_PVTR_File(DirName, FileName; field=nothing) +function Read_LaMEM_PVTR_File(DirName::String, FileName::String; field=nothing) CurDir = pwd(); # change to directory @@ -196,9 +196,9 @@ end """ - FileNames, Time = readPVD(FileName::String) + FileNames, Time, Timestep = readPVD(FileName::String) -This reads a PVD file & returns the timesteps and corresponding filenames +This reads a PVD file & returns the `FileNames`, `Time` and `Timesteps` """ function readPVD(FileName::String) @@ -206,18 +206,25 @@ function readPVD(FileName::String) start_line = findall(lines .== "")[1] + 1 end_line = findall(lines .== "")[1] - 1 - FileNames = []; - Time = []; + FileNames = Vector{String}() + Time = Vector{Float64}() + Timestep = Vector{Int64}() for i=start_line:end_line line = lines[i] time = split(line)[2]; time = parse(Float64,time[11:end-1]) file = split(line)[3]; file = String.(file[7:end-3]); - FileNames = [FileNames; file] - Time = [Time; time] + FileNames = push!(FileNames, file) + Time = push!(Time, time) + + # retrieve the timestep + file_name = split(file,Base.Filesystem.path_separator)[1]; + timestep = parse(Int64,split(file_name,"_")[2]); + Timestep = push!(Timestep, timestep) + end - return FileNames, Time + return FileNames, Time, Timestep end @@ -322,4 +329,49 @@ function Read_LaMEM_PVTS_File(DirName, FileName; field=nothing) data_output = CartData(X,Y,Z, data_fields) return data_output +end + + +""" + +This reads a LaMEM timestep. + +Input Arguments: +- `FileName`: name of the simulation, w/out extension +- `DirName`: name of the main directory (i.e. where the `*.pvd` files are located) + +""" +function Read_LaMEM_timestep(FileName::String, TimeStep::Int64=0, DirName::String=""; field=nothing, phase=false, surf=false) + + FileNames, Time, Timestep = Read_LaMEM_simulation(FileName, DirName; phase=phase, surf=surf); + @show Timestep + ind = findall(Timestep.==TimeStep) + + if isempty(ind) + error("this timestep does not exist") + end + + + +end + + +""" + FileNames, Time, Timestep = Read_LaMEM_simulation(FileName::String, DirName::String=""; phase=false, surf=false) + +Reads a LaMEM simulation `FileName` in directory `DirName` and returns the timesteps, times and filenames of that simulation. +""" +function Read_LaMEM_simulation(FileName::String, DirName::String=""; phase=false, surf=false) + + if phase==true + pvd_file=FileName*"_phase.pvd" + elseif surf==true + pvd_file=FileName*"_surf.pvd" + else + pvd_file=FileName*".pvd" + end + + FileNames, Time, Timestep = readPVD(joinpath(DirName,pvd_file)) + + return FileNames, Time, Timestep end \ No newline at end of file diff --git a/src/run_lamem_save_grid.jl b/src/run_lamem_save_grid.jl index 5b2555c6..c17673f0 100644 --- a/src/run_lamem_save_grid.jl +++ b/src/run_lamem_save_grid.jl @@ -4,23 +4,6 @@ # Note: This downloads the BinaryBuilder version of LaMEM, which is not necessarily the latest version of LaMEM # (or the same as the current repository), since we have to manually update the builds. -""" - run_lamem_save_grid(ParamFile::String, cores::Int64=1) -This calls LaMEM simulation, for using the parameter file `ParamFile` -and creates processor paritioning file "ProcessorPartitioning_`cores`cpu_X.Y.Z.bin" for `cores` number of cores. -# Example: -The first step is to ensure that `LaMEM_jll` is installed on your system. You only need to do this once, or once LaMEM_jll is updated. -```julia -julia> import Pkg -julia> Pkg.add("LaMEM_jll") -``` -Next you can call LaMEM with: -```julia -julia> ParamFile="../../input_models/BuildInSetups/FallingBlock_Multigrid.dat"; -julia> run_lamem_save_grid(ParamFile, 2) -``` -""" - function run_lamem_with_log(ParamFile::String, cores::Int64=1, args::String="") @@ -48,7 +31,7 @@ function JuliaStringToArray(input) return arr end -function get_line_contatinig(stringarray::Vector{SubString{String}}, lookfor::String) +function get_line_containing(stringarray::Vector{SubString{String}}, lookfor::String) for line in stringarray @@ -59,12 +42,31 @@ function get_line_contatinig(stringarray::Vector{SubString{String}}, lookfor::St end end +""" + run_lamem_save_grid(ParamFile::String, cores::Int64=1) +This calls LaMEM simulation, for using the parameter file `ParamFile` +and creates processor paritioning file "ProcessorPartitioning_`cores`cpu_X.Y.Z.bin" for `cores` number of cores. +# Example: +The first step is to ensure that `LaMEM_jll` is installed on your system. You only need to do this once, or once LaMEM_jll is updated. +```julia +julia> import Pkg +julia> Pkg.add("LaMEM_jll") +``` +Next you can call LaMEM with: +```julia +julia> ParamFile="../../input_models/BuildInSetups/FallingBlock_Multigrid.dat"; +julia> run_lamem_save_grid(ParamFile, 2) +``` +""" function run_lamem_save_grid(ParamFile::String, cores::Int64=1) - if cores==1 return print("No partitioning file required for 1 core model setup \n") end + if cores==1 + return print("No partitioning file required for 1 core model setup \n") + end + ParamFile = abspath(ParamFile) logoutput = run_lamem_with_log(ParamFile, cores,"-mode save_grid" ) arr = JuliaStringToArray(logoutput) - foundline = get_line_contatinig(arr,"Processor grid [nx, ny, nz] : ") + foundline = get_line_containing(arr,"Processor grid [nx, ny, nz] : ") foundline = join(map(x -> isspace(foundline[x]) ? "" : foundline[x], 1:length(foundline))) sprtlftbrkt = split(foundline,"[") sprtrghtbrkt = split(sprtlftbrkt[3],"]") diff --git a/src/utils_IO.jl b/src/utils_IO.jl new file mode 100644 index 00000000..40922fe7 --- /dev/null +++ b/src/utils_IO.jl @@ -0,0 +1,80 @@ +using Glob + +export clean_directory, changefolder + +""" + clean_directory(DirName) + +Removes all LaMEM timesteps & `*.pvd` files from the directory `DirName` + +""" +function clean_directory(DirName="./") + + CurDir = pwd(); + + # change to directory + cd(DirName) + + # pvd files + for f in glob("*.pvd") + rm(f) + end + + # vts files + for f in glob("*.vts") + rm(f) + end + + #timestep directories + for f in glob("Timestep*") + rm(f, recursive=true, force=true) + end + + + cd(CurDir) +end + + +""" + changefolder() + +Starts a GUI on Windowss or Mac machines, which allows you to change our working directory +""" +function changefolder() + if Sys.iswindows() + command = """ + Function Get-Folder(\$initialDirectory) { + [System.Reflection.Assembly]::LoadWithPartialName("System.windows.forms")|Out-Null + + \$foldername = New-Object System.Windows.Forms.FolderBrowserDialog + \$foldername.Description = "Select a folder" + \$foldername.rootfolder = "MyComputer" + + if(\$foldername.ShowDialog() -eq "OK") + { + \$folder += \$foldername.SelectedPath + } + return \$folder + } + + Get-Folder + """ + cd(chomp(read(`powershell -Command $command`, String))) + println(pwd()) + elseif Sys.isapple() + command = """ + try + set af to (choose folder with prompt "Folder?") + set result to POSIX path of af + on error + beep + set result to "$(pwd())" + end + result + """ + cd(chomp(read(`osascript -e $command`, String))) + println(pwd()) + else + exit() + end +end diff --git a/src/utils.jl b/src/utils_Run.jl similarity index 53% rename from src/utils.jl rename to src/utils_Run.jl index 62b3c7f7..847382d7 100644 --- a/src/utils.jl +++ b/src/utils_Run.jl @@ -1,8 +1,4 @@ -using LaMEM_jll -using Glob - -export remove_popup_messages_mac, show_paths_LaMEM, clean_directory, changefolder - +export remove_popup_messages_mac, show_paths_LaMEM """ remove_popup_messages_mac() @@ -74,81 +70,3 @@ function show_paths_LaMEM() return nothing end - -""" - clean_directory(DirName) - -Removes all LaMEM timesteps & `*.pvd` files from the directory `DirName` - -""" -function clean_directory(DirName="./") - - CurDir = pwd(); - - # change to directory - cd(DirName) - - # pvd files - for f in glob("*.pvd") - rm(f) - end - - # vts files - for f in glob("*.vts") - rm(f) - end - - #timestep directories - for f in glob("Timestep*") - rm(f, recursive=true, force=true) - end - - - cd(CurDir) - -end - - -""" - changefolder() - -Starts a GUI on Windowss or Mac machines, which allows you to change our working directory -""" -function changefolder() - if Sys.iswindows() - command = """ - Function Get-Folder(\$initialDirectory) { - [System.Reflection.Assembly]::LoadWithPartialName("System.windows.forms")|Out-Null - - \$foldername = New-Object System.Windows.Forms.FolderBrowserDialog - \$foldername.Description = "Select a folder" - \$foldername.rootfolder = "MyComputer" - - if(\$foldername.ShowDialog() -eq "OK") - { - \$folder += \$foldername.SelectedPath - } - return \$folder - } - - Get-Folder - """ - cd(chomp(read(`powershell -Command $command`, String))) - println(pwd()) - elseif Sys.isapple() - command = """ - try - set af to (choose folder with prompt "Folder?") - set result to POSIX path of af - on error - beep - set result to "$(pwd())" - end - result - """ - cd(chomp(read(`osascript -e $command`, String))) - println(pwd()) - else - exit() - end -end From 4db4054fd6af36ff00f68858a54b91e2155d4efa Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Sat, 8 Apr 2023 22:45:18 +0200 Subject: [PATCH 09/18] this simplifies reading files & figuring out which data is stored in the LaMEM output files --- src/LaMEM.jl | 4 +- src/read_timestep.jl | 149 +++++++++++++++++++++++++++++++++---------- 2 files changed, 117 insertions(+), 36 deletions(-) diff --git a/src/LaMEM.jl b/src/LaMEM.jl index 328bd47d..a5c5a704 100644 --- a/src/LaMEM.jl +++ b/src/LaMEM.jl @@ -3,8 +3,8 @@ module LaMEM # Functions to read LaMEM output include("IO.jl") using .IO -export Read_LaMEM_PVTR_File, Read_LaMEM_PVTS_File, field_names, readPVD, Read_LaMEM_PVTU_File -export Read_LaMEM_simulation, Read_LaMEM_timestep +export Read_LaMEM_PVTR_File, Read_LaMEM_PVTS_File, Read_LaMEM_PVTU_File +export Read_LaMEM_simulation, Read_LaMEM_timestep, Read_LaMEM_fieldnames # Functions include("Run.jl") diff --git a/src/read_timestep.jl b/src/read_timestep.jl index bb8fab66..eedbe405 100644 --- a/src/read_timestep.jl +++ b/src/read_timestep.jl @@ -6,8 +6,8 @@ using GeophysicalModelGenerator: CartData, XYZGrid using Glob, ReadVTK -export Read_LaMEM_PVTR_File, Read_LaMEM_PVTS_File, field_names, Read_LaMEM_PVTU_File -export Read_LaMEM_simulation, Read_LaMEM_timestep +export Read_LaMEM_PVTR_File, Read_LaMEM_PVTS_File, Read_LaMEM_PVTU_File +export Read_LaMEM_simulation, Read_LaMEM_timestep, Read_LaMEM_fieldnames """ output, isCell = ReadField_3D_pVTR(data, FieldName::String) @@ -125,10 +125,19 @@ function ReadField_3D_pVTU(pvtu, FieldName) return data_out end + # The FileName can contain a directory as well; deal with that here +function split_path_name(DirName_base::String, FileName::String) + FullName = joinpath(DirName_base,FileName) + id = findlast(Base.Filesystem.path_separator, FullName)[1]; + DirName = FullName[1:id-1] + File = FullName[id+1:end] + + return DirName, File +end """ - data_output = Read_LaMEM_PVTR_File(DirName, FileName; field=nothing) + data_output = Read_LaMEM_PVTR_File(DirName, FileName; fields=nothing) Reads a 3D LaMEM timestep from VTR file `FileName`, located in directory `DirName`. By default, it will read all fields. If you want you can only read a specific `field`. See the function `fieldnames` to get a list with all available fields in the file. @@ -136,16 +145,16 @@ By default, it will read all fields. If you want you can only read a specific `f It will return `data_output` which is a `CartData` output structure. """ -function Read_LaMEM_PVTR_File(DirName::String, FileName::String; field=nothing) +function Read_LaMEM_PVTR_File(DirName_base::String, FileName::String; fields=nothing) CurDir = pwd(); + DirName, File = split_path_name(DirName_base, FileName) + # change to directory cd(DirName) # read data from parallel rectilinear grid - pvtk = PVTKFile(FileName) - - cd(CurDir) + pvtk = PVTKFile(File) # Fields stored in this data file: names = keys(get_point_data(pvtk)) @@ -154,7 +163,7 @@ function Read_LaMEM_PVTR_File(DirName::String, FileName::String; field=nothing) end isCell = true - if isnothing(field) + if isnothing(fields) # read all data in the file data_fields = NamedTuple(); @@ -164,15 +173,20 @@ function Read_LaMEM_PVTR_File(DirName::String, FileName::String; field=nothing) end else - # read just a single data set - ind = findall(contains.(field, names)) - # check that it exists - if isempty(ind) - error("the field $field does not exist in the data file") - else - data_fields, isCell = ReadField_3D_pVTR(data, field) + # read the data sets specified + data_fields = NamedTuple(); + for field in fields + ind = findall(contains.(field, names)) + # check that it exists + if isempty(ind) + error("the field $field does not exist in the data file") + else + dat, isCell = ReadField_3D_pVTR(pvtk, field) + end + data_fields = merge(data_fields,dat) end end + cd(CurDir) coords_read = get_coordinates(pvtk) @@ -255,6 +269,30 @@ function Read_LaMEM_PVTU_File(DirName, FileName; field=nothing) end isCell = true + #= + if isnothing(fields) + # read all data in the file + data_fields = NamedTuple(); + for FieldName in names + dat, isCell = ReadField_3D_pVTR(pvtk, FieldName); + data_fields = merge(data_fields,dat) + end + + else + # read the data sets specified + data_fields = NamedTuple(); + for field in fields + ind = findall(contains.(field, names)) + # check that it exists + if isempty(ind) + error("the field $field does not exist in the data file") + else + dat, isCell = ReadField_3D_pVTR(pvtk, field) + end + data_fields = merge(data_fields,dat) + end + end + =# if isnothing(field) # read all data in the file data_fields = NamedTuple(); @@ -295,22 +333,20 @@ By default, it will read all fields. If you want you can only read a specific `f It will return `data_output` which is a `CartData` output structure. """ -function Read_LaMEM_PVTS_File(DirName, FileName; field=nothing) +function Read_LaMEM_PVTS_File(DirName_base::String, FileName::String; fields=nothing) CurDir = pwd(); + DirName, File = split_path_name(DirName_base, FileName) + # read data from parallel rectilinear grid cd(DirName) - pvts = PVTKFile(FileName) + pvts = PVTKFile(File) cd(CurDir) # Fields stored in this data file: - name = keys(get_point_data(pvts)) - if !isnothing(field) - if any(occursin.(name,field)) - name = tuple(field); - else - error("the field $field does not exist in the data file which has fields: $(name)") - end + names = keys(get_point_data(pvts)) + if isempty(names) + names = keys(get_cell_data(pvts)) end # Read coordinates @@ -318,8 +354,8 @@ function Read_LaMEM_PVTS_File(DirName, FileName; field=nothing) # read all data in the file data_fields = NamedTuple(); - for FieldName in name - dat, _ = LaMEM.ReadField_3D_pVTS(pvts, FieldName); + for FieldName in names + dat, _ = ReadField_3D_pVTS(pvts, FieldName); if length(dat[1])>length(X) dat_tup = ntuple(i->dat[1][i,:,:,:], size(dat[1],1)) dat = NamedTuple{keys(dat)}(tuple(dat_tup)) @@ -333,26 +369,41 @@ end """ + data, time = Read_LaMEM_timestep(FileName::String, TimeStep::Int64=0, DirName::String=""; fields=nothing, phase=false, surf=false, last=false) This reads a LaMEM timestep. Input Arguments: - `FileName`: name of the simulation, w/out extension +- `Timestep`: timestep to be read, unless `last=true` in which case we read the last one - `DirName`: name of the main directory (i.e. where the `*.pvd` files are located) +- `fields`: Tuple with optional fields; if not specified all will be loaded +- `phase`: Loads the phase information of LaMEM if true +- `surf`: Loads the free surface of LaMEM if true +- `last`: Loads the last timestep + +Output: +- `data`: Cartesian data struct with LaMEM output +- `time`: The time of the timestep """ -function Read_LaMEM_timestep(FileName::String, TimeStep::Int64=0, DirName::String=""; field=nothing, phase=false, surf=false) +function Read_LaMEM_timestep(FileName::String, TimeStep::Int64=0, DirName::String=""; fields=nothing, phase=false, surf=false, last=false) - FileNames, Time, Timestep = Read_LaMEM_simulation(FileName, DirName; phase=phase, surf=surf); - @show Timestep + Timestep, FileNames, Time = Read_LaMEM_simulation(FileName, DirName; phase=phase, surf=surf); + ind = findall(Timestep.==TimeStep) - - if isempty(ind) - error("this timestep does not exist") - end + if last==true; ind = length(Time); end + if isempty(ind); error("this timestep does not exist"); end + # Read file + if surf==true + data = Read_LaMEM_PVTS_File(DirName, FileNames[ind[1]], fields=fields) + else + data = Read_LaMEM_PVTR_File(DirName, FileNames[ind[1]], fields=fields) + end + return data, Time[ind] end @@ -373,5 +424,35 @@ function Read_LaMEM_simulation(FileName::String, DirName::String=""; phase=false FileNames, Time, Timestep = readPVD(joinpath(DirName,pvd_file)) - return FileNames, Time, Timestep + return Timestep, FileNames, Time +end + +""" + Read_LaMEM_fieldnames(FileName::String, DirName_base::String=""; phase=false, surf=false) + +Returns the names of the datasets stored in `FileName` +""" +function Read_LaMEM_fieldnames(FileName::String, DirName_base::String=""; phase=false, surf=false) + + _, FileNames, _ = Read_LaMEM_simulation(FileName, DirName_base; phase=phase, surf=surf); + + # Read file + DirName, File = split_path_name(DirName_base, FileNames[1]) + + # change to directory + cur_dir = pwd(); + + # read data from parallel rectilinear grid + cd(DirName) + pvtk = PVTKFile(File) + cd(cur_dir) + + # Fields stored in this data file: + if phase==false + names = keys(get_point_data(pvtk)) + else + names = keys(get_cell_data(pvtk)) + end + + return names end \ No newline at end of file From 69c7d8f54240a538adc667fc2a9a1403b0ffc50e Mon Sep 17 00:00:00 2001 From: bkaus Date: Sat, 8 Apr 2023 23:16:44 +0200 Subject: [PATCH 10/18] cleanup --- src/LaMEM.jl | 23 ----------------------- 1 file changed, 23 deletions(-) diff --git a/src/LaMEM.jl b/src/LaMEM.jl index a5c5a704..486f7eb9 100644 --- a/src/LaMEM.jl +++ b/src/LaMEM.jl @@ -12,27 +12,4 @@ using .Run export run_lamem, run_lamem_save_grid, mpiexec export remove_popup_messages_mac, show_paths_LaMEM - - -#= -using LaMEM_jll - -# load the correct MPI -const mpiexec = if isdefined(LaMEM_jll,:MPICH_jll) - LaMEM_jll.MPICH_jll.mpiexec() -elseif isdefined(LaMEM_jll,:MicrosoftMPI_jll) - LaMEM_jll.MicrosoftMPI_jll.mpiexec() -else - nothing -end - - - -=# - - -#include("read_timestep.jl") -#include("utils.jl") - - end # module From e538121acb27afc521f9e49de909260ef863c6c6 Mon Sep 17 00:00:00 2001 From: bkaus Date: Sat, 8 Apr 2023 23:54:34 +0200 Subject: [PATCH 11/18] make tests work --- Project.toml | 3 ++- src/LaMEM.jl | 1 + test/read_timestep.jl | 29 +++++++++++++++-------------- test/runLaMEM.jl | 13 ++++++------- 4 files changed, 24 insertions(+), 22 deletions(-) diff --git a/Project.toml b/Project.toml index 8aec278b..1bed6337 100644 --- a/Project.toml +++ b/Project.toml @@ -12,8 +12,9 @@ ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3" [compat] GeophysicalModelGenerator = "0.4" Glob = "1" -LaMEM_jll = "1.2" +LaMEM_jll = "1.2.3" julia = "1.6" +ReadVTK = "0.1.5" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/src/LaMEM.jl b/src/LaMEM.jl index 486f7eb9..ed0e177b 100644 --- a/src/LaMEM.jl +++ b/src/LaMEM.jl @@ -5,6 +5,7 @@ include("IO.jl") using .IO export Read_LaMEM_PVTR_File, Read_LaMEM_PVTS_File, Read_LaMEM_PVTU_File export Read_LaMEM_simulation, Read_LaMEM_timestep, Read_LaMEM_fieldnames +export clean_directory, changefolder # Functions include("Run.jl") diff --git a/test/read_timestep.jl b/test/read_timestep.jl index 4e3f4de7..0b24a35b 100644 --- a/test/read_timestep.jl +++ b/test/read_timestep.jl @@ -4,36 +4,37 @@ using LaMEM @testset "read LaMEM output" begin # Read a timestep - FileName="FB_multigrid.pvtr" - DirName = "Timestep_00000001_6.72970343e+00" - - data = Read_LaMEM_PVTR_File(DirName, FileName) + FileName="FB_multigrid" + Timestep = 1 + data, time = Read_LaMEM_timestep(FileName,Timestep) @test sum(data.fields.phase) ≈ 736.36414f0 @test sum(data.fields.strain_rate[1,:,:,:]) ≈ -0.019376338f0 - fields = field_names(DirName, FileName) - + fields = Read_LaMEM_fieldnames(FileName) + # with cell-data - FileName="FB_multigrid_phase.pvtr" + FileName="FB_multigrid" DirName = "Timestep_00000001_6.72970343e+00" - data = Read_LaMEM_PVTR_File(DirName, FileName) + Timestep = 1 + data, time = Read_LaMEM_timestep(FileName,Timestep, phase=true) + @test sum(data.fields.phase) == 19822 # read subduction setup - data = Read_LaMEM_PVTR_File("Timestep_00000001_5.50000000e-02", "Subduction2D_FreeSlip_direct.pvtr") + data, time = Read_LaMEM_timestep("Subduction2D_FreeSlip_direct",1) @test sum(data.fields.density) ≈ 1.60555f8 # Read PVD file - FileNames, Time = readPVD("Subduction2D_FreeSlip_direct.pvd") + Timestep, FileNames, Time = Read_LaMEM_simulation("Subduction2D_FreeSlip_direct") @test Time[2] ≈ 0.055 # Read passive tracers - data = Read_LaMEM_PVTU_File("Timestep_00000010_1.09635548e+00", "PlumeLithosphereInteraction_passive_tracers.pvtu") - @test data.z[100] ≈ -298.4531f0 - @test data.fields.Temperature[100] ≈ 1350.0f0 + #data = Read_LaMEM_PVTU_File("Timestep_00000010_1.09635548e+00", "PlumeLithosphereInteraction_passive_tracers.pvtu") + #@test data.z[100] ≈ -298.4531f0 + #@test data.fields.Temperature[100] ≈ 1350.0f0 # Read surface data - data = Read_LaMEM_PVTS_File("Timestep_00000005_7.59607376e-02", "Subduction2D_FreeSurface_direct_surf.pvts") + data, time = Read_LaMEM_timestep("Subduction2D_FreeSurface_direct",5, surf=true) @test data.z[100] ≈ 0.6830405f0 @test sum(data.fields.topography[:,1,1]) ≈ 1.2634416f0 end diff --git a/test/runLaMEM.jl b/test/runLaMEM.jl index fb7e8e00..9f43dd40 100644 --- a/test/runLaMEM.jl +++ b/test/runLaMEM.jl @@ -29,18 +29,17 @@ using LaMEM out = run_lamem(ParamFile, 1, "-nstep_max 2") # 1 core @test isnothing(out) + if !Sys.iswindows() + out = run_lamem(ParamFile, 2, "-nstep_max 5") # 2 cores (mumps) + @test isnothing(out) + end + # Try free surface ParamFile="input_files/Subduction2D_FreeSurface_DirectSolver.dat"; out = run_lamem(ParamFile, 4, "-nstep_max 5") # 4 core @test isnothing(out) - - if !Sys.iswindows() - out = run_lamem(ParamFile, 2, "-nstep_max 2") # 2 cores (mumps) - @test isnothing(out) - end - - if !Sys.isapple() + if !Sys.isapple() && 1==0 # Note: superlu_dist uses a combination of openMP parallelization on a node and MPI between nodes. # If you have a server with a lot of cores AND run this with >1 core, this may clash # In that case it is better to run it with 1 MPI task but set the environmental variables below accordingly From 416d4a1ea34ebcffe7cc81284bf74964fec8f044 Mon Sep 17 00:00:00 2001 From: bkaus Date: Sun, 9 Apr 2023 00:22:20 +0200 Subject: [PATCH 12/18] make it work for passive tracers --- src/read_timestep.jl | 55 ++++++++++++++++++------------------------- test/read_timestep.jl | 6 ++--- 2 files changed, 26 insertions(+), 35 deletions(-) diff --git a/src/read_timestep.jl b/src/read_timestep.jl index eedbe405..921af98c 100644 --- a/src/read_timestep.jl +++ b/src/read_timestep.jl @@ -243,7 +243,7 @@ end """ - data_output = Read_LaMEM_PVTU_File(DirName, FileName; field=nothing) + data_output = Read_LaMEM_PVTU_File(DirName, FileName; fields=nothing) Reads a 3D LaMEM timestep from VTU file `FileName`, located in directory `DirName`. Typically this is done to read passive tracers back into julia. By default, it will read all fields. If you want you can only read a specific `field`. See the function `fieldnames` to get a list with all available fields in the file. @@ -251,14 +251,16 @@ By default, it will read all fields. If you want you can only read a specific `f It will return `data_output` which is a `CartData` output structure. """ -function Read_LaMEM_PVTU_File(DirName, FileName; field=nothing) +function Read_LaMEM_PVTU_File(DirName_base, FileName; fields=nothing) CurDir = pwd(); + DirName, File = split_path_name(DirName_base, FileName) + # change to directory cd(DirName) # read data from parallel rectilinear grid - pvtu = PVTKFile(FileName) + pvtu = PVTKFile(File) cd(CurDir) @@ -269,12 +271,11 @@ function Read_LaMEM_PVTU_File(DirName, FileName; field=nothing) end isCell = true - #= if isnothing(fields) # read all data in the file data_fields = NamedTuple(); for FieldName in names - dat, isCell = ReadField_3D_pVTR(pvtk, FieldName); + dat = ReadField_3D_pVTU(pvtu, FieldName); data_fields = merge(data_fields,dat) end @@ -287,30 +288,11 @@ function Read_LaMEM_PVTU_File(DirName, FileName; field=nothing) if isempty(ind) error("the field $field does not exist in the data file") else - dat, isCell = ReadField_3D_pVTR(pvtk, field) + dat, isCell = ReadField_3D_pVTU(pvtk, field) end data_fields = merge(data_fields,dat) end end - =# - if isnothing(field) - # read all data in the file - data_fields = NamedTuple(); - for FieldName in names - dat = ReadField_3D_pVTU(pvtu, FieldName); - data_fields = merge(data_fields,dat) - end - - else - # read just a single data set - ind = findall(contains.(field, names)) - # check that it exists - if isempty(ind) - error("the field $field does not exist in the data file") - else - data_fields = ReadField_3D_pVTU(data, field) - end - end points = get_points(pvtu) @@ -380,6 +362,7 @@ Input Arguments: - `fields`: Tuple with optional fields; if not specified all will be loaded - `phase`: Loads the phase information of LaMEM if true - `surf`: Loads the free surface of LaMEM if true +- `passive_tracers`: Loads passive tracers if true - `last`: Loads the last timestep Output: @@ -387,9 +370,9 @@ Output: - `time`: The time of the timestep """ -function Read_LaMEM_timestep(FileName::String, TimeStep::Int64=0, DirName::String=""; fields=nothing, phase=false, surf=false, last=false) +function Read_LaMEM_timestep(FileName::String, TimeStep::Int64=0, DirName::String=pwd(); fields=nothing, phase=false, surf=false, passive_tracers=false, last=false) - Timestep, FileNames, Time = Read_LaMEM_simulation(FileName, DirName; phase=phase, surf=surf); + Timestep, FileNames, Time = Read_LaMEM_simulation(FileName, DirName; phase=phase, surf=surf, passive_tracers=passive_tracers); ind = findall(Timestep.==TimeStep) @@ -399,6 +382,8 @@ function Read_LaMEM_timestep(FileName::String, TimeStep::Int64=0, DirName::Strin # Read file if surf==true data = Read_LaMEM_PVTS_File(DirName, FileNames[ind[1]], fields=fields) + elseif passive_tracers==true + data = Read_LaMEM_PVTU_File(DirName, FileNames[ind[1]], fields=fields) else data = Read_LaMEM_PVTR_File(DirName, FileNames[ind[1]], fields=fields) end @@ -408,16 +393,18 @@ end """ - FileNames, Time, Timestep = Read_LaMEM_simulation(FileName::String, DirName::String=""; phase=false, surf=false) + FileNames, Time, Timestep = Read_LaMEM_simulation(FileName::String, DirName::String=""; phase=false, surf=false, passive_tracers=false) Reads a LaMEM simulation `FileName` in directory `DirName` and returns the timesteps, times and filenames of that simulation. """ -function Read_LaMEM_simulation(FileName::String, DirName::String=""; phase=false, surf=false) +function Read_LaMEM_simulation(FileName::String, DirName::String=""; phase=false, surf=false, passive_tracers=false) if phase==true pvd_file=FileName*"_phase.pvd" elseif surf==true pvd_file=FileName*"_surf.pvd" + elseif passive_tracers==true + pvd_file=FileName*"_passive_tracers.pvd" else pvd_file=FileName*".pvd" end @@ -428,11 +415,11 @@ function Read_LaMEM_simulation(FileName::String, DirName::String=""; phase=false end """ - Read_LaMEM_fieldnames(FileName::String, DirName_base::String=""; phase=false, surf=false) + Read_LaMEM_fieldnames(FileName::String, DirName_base::String=""; phase=false, surf=false, tracers=false) Returns the names of the datasets stored in `FileName` """ -function Read_LaMEM_fieldnames(FileName::String, DirName_base::String=""; phase=false, surf=false) +function Read_LaMEM_fieldnames(FileName::String, DirName_base::String=""; phase=false, surf=false, tracers=false) _, FileNames, _ = Read_LaMEM_simulation(FileName, DirName_base; phase=phase, surf=surf); @@ -444,7 +431,11 @@ function Read_LaMEM_fieldnames(FileName::String, DirName_base::String=""; phase= # read data from parallel rectilinear grid cd(DirName) - pvtk = PVTKFile(File) + if !tracers + pvtk = PVTKFile(File) + else + pvtk = PVTUFile(File) + end cd(cur_dir) # Fields stored in this data file: diff --git a/test/read_timestep.jl b/test/read_timestep.jl index 0b24a35b..e56feb5f 100644 --- a/test/read_timestep.jl +++ b/test/read_timestep.jl @@ -29,9 +29,9 @@ using LaMEM @test Time[2] ≈ 0.055 # Read passive tracers - #data = Read_LaMEM_PVTU_File("Timestep_00000010_1.09635548e+00", "PlumeLithosphereInteraction_passive_tracers.pvtu") - #@test data.z[100] ≈ -298.4531f0 - #@test data.fields.Temperature[100] ≈ 1350.0f0 + data, time = Read_LaMEM_timestep("PlumeLithosphereInteraction",10, passive_tracers=true) + @test data.z[100] ≈ -298.4531f0 + @test data.fields.Temperature[100] ≈ 1350.0f0 # Read surface data data, time = Read_LaMEM_timestep("Subduction2D_FreeSurface_direct",5, surf=true) From e9ad1bb6fb192326daa4e44ad9c06e4277f0b750 Mon Sep 17 00:00:00 2001 From: bkaus Date: Sun, 9 Apr 2023 13:19:05 +0200 Subject: [PATCH 13/18] test on apple --- .../Subduction2D_FreeSurface_DirectSolver.dat | 2 +- test/read_timestep.jl | 10 ++++++---- test/runLaMEM.jl | 10 ++++++---- 3 files changed, 13 insertions(+), 9 deletions(-) diff --git a/test/input_files/Subduction2D_FreeSurface_DirectSolver.dat b/test/input_files/Subduction2D_FreeSurface_DirectSolver.dat index 6b1d1d07..2728a177 100644 --- a/test/input_files/Subduction2D_FreeSurface_DirectSolver.dat +++ b/test/input_files/Subduction2D_FreeSurface_DirectSolver.dat @@ -78,7 +78,7 @@ # Solver options #=============================================================================== SolverType = direct # solver [direct or multigrid] - DirectSolver = mumps # mumps/superlu_dist/pastix + #DirectSolver = mumps # mumps/superlu_dist/pastix DirectPenalty = 1e5 diff --git a/test/read_timestep.jl b/test/read_timestep.jl index e56feb5f..ae4891cd 100644 --- a/test/read_timestep.jl +++ b/test/read_timestep.jl @@ -33,9 +33,11 @@ using LaMEM @test data.z[100] ≈ -298.4531f0 @test data.fields.Temperature[100] ≈ 1350.0f0 - # Read surface data - data, time = Read_LaMEM_timestep("Subduction2D_FreeSurface_direct",5, surf=true) - @test data.z[100] ≈ 0.6830405f0 - @test sum(data.fields.topography[:,1,1]) ≈ 1.2634416f0 + if !Sys.isapple() # broken on mac for LaMEM_jll 1.2.3 (should be fixed in next release) + # Read surface data + data, time = Read_LaMEM_timestep("Subduction2D_FreeSurface_direct",5, surf=true) + @test data.z[100] ≈ 0.6830405f0 + @test sum(data.fields.topography[:,1,1]) ≈ 1.2634416f0 + end end diff --git a/test/runLaMEM.jl b/test/runLaMEM.jl index 9f43dd40..fa726991 100644 --- a/test/runLaMEM.jl +++ b/test/runLaMEM.jl @@ -34,10 +34,12 @@ using LaMEM @test isnothing(out) end - # Try free surface - ParamFile="input_files/Subduction2D_FreeSurface_DirectSolver.dat"; - out = run_lamem(ParamFile, 4, "-nstep_max 5") # 4 core - @test isnothing(out) + if !Sys.isapple() # broken on mac for LaMEM_jll 1.2.3 (should be fixed in next release) + # Try free surface + ParamFile="input_files/Subduction2D_FreeSurface_DirectSolver.dat"; + out = run_lamem(ParamFile, 1, "-nstep_max 5") # 4 core + @test isnothing(out) + end if !Sys.isapple() && 1==0 # Note: superlu_dist uses a combination of openMP parallelization on a node and MPI between nodes. From ca3b5f7aa7ef6e3e78f322bdbe4a0cc95c01fc61 Mon Sep 17 00:00:00 2001 From: bkaus Date: Sun, 9 Apr 2023 13:28:55 +0200 Subject: [PATCH 14/18] update --- test/read_timestep.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/read_timestep.jl b/test/read_timestep.jl index ae4891cd..b36806a5 100644 --- a/test/read_timestep.jl +++ b/test/read_timestep.jl @@ -36,8 +36,8 @@ using LaMEM if !Sys.isapple() # broken on mac for LaMEM_jll 1.2.3 (should be fixed in next release) # Read surface data data, time = Read_LaMEM_timestep("Subduction2D_FreeSurface_direct",5, surf=true) - @test data.z[100] ≈ 0.6830405f0 - @test sum(data.fields.topography[:,1,1]) ≈ 1.2634416f0 + @test data.z[100] ≈ 0.68236357f0 + @test sum(data.fields.topography[:,1,1]) ≈ 1.2645866f0 end end From 396f5a241ba34137b1d62c095e14864764787738 Mon Sep 17 00:00:00 2001 From: bkaus Date: Sun, 9 Apr 2023 16:31:03 +0200 Subject: [PATCH 15/18] debug for windows --- test/read_timestep.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/read_timestep.jl b/test/read_timestep.jl index b36806a5..9129625b 100644 --- a/test/read_timestep.jl +++ b/test/read_timestep.jl @@ -6,6 +6,7 @@ using LaMEM # Read a timestep FileName="FB_multigrid" Timestep = 1 + @show readdir() data, time = Read_LaMEM_timestep(FileName,Timestep) @test sum(data.fields.phase) ≈ 736.36414f0 @test sum(data.fields.strain_rate[1,:,:,:]) ≈ -0.019376338f0 From 349fd5ac1c242f97862762e0e177392da25b52d9 Mon Sep 17 00:00:00 2001 From: bkaus Date: Mon, 10 Apr 2023 10:37:25 +0200 Subject: [PATCH 16/18] change file seperator (also on win) --- src/read_timestep.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/read_timestep.jl b/src/read_timestep.jl index 921af98c..51b54235 100644 --- a/src/read_timestep.jl +++ b/src/read_timestep.jl @@ -232,7 +232,9 @@ function readPVD(FileName::String) Time = push!(Time, time) # retrieve the timestep - file_name = split(file,Base.Filesystem.path_separator)[1]; + #file_name = split(file,Base.Filesystem.path_separator)[1]; + file_name = split(file,"/")[1]; + timestep = parse(Int64,split(file_name,"_")[2]); Timestep = push!(Timestep, timestep) From d5b1bf66c092931e1955a8d83877eff608b57d74 Mon Sep 17 00:00:00 2001 From: bkaus Date: Mon, 10 Apr 2023 10:44:23 +0200 Subject: [PATCH 17/18] next attempt --- src/read_timestep.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/read_timestep.jl b/src/read_timestep.jl index 51b54235..ce3bf672 100644 --- a/src/read_timestep.jl +++ b/src/read_timestep.jl @@ -128,7 +128,9 @@ end # The FileName can contain a directory as well; deal with that here function split_path_name(DirName_base::String, FileName::String) FullName = joinpath(DirName_base,FileName) - id = findlast(Base.Filesystem.path_separator, FullName)[1]; + #id = findlast(Base.Filesystem.path_separator, FullName)[1]; + id = findlast("/", FullName)[1]; + DirName = FullName[1:id-1] File = FullName[id+1:end] @@ -148,7 +150,9 @@ It will return `data_output` which is a `CartData` output structure. function Read_LaMEM_PVTR_File(DirName_base::String, FileName::String; fields=nothing) CurDir = pwd(); + @show DirName_base FileName DirName, File = split_path_name(DirName_base, FileName) + @show DirName File # change to directory cd(DirName) @@ -232,8 +236,8 @@ function readPVD(FileName::String) Time = push!(Time, time) # retrieve the timestep - #file_name = split(file,Base.Filesystem.path_separator)[1]; - file_name = split(file,"/")[1]; + file_name = split(file,Base.Filesystem.path_separator)[1]; + #file_name = split(file,"/")[1]; timestep = parse(Int64,split(file_name,"_")[2]); Timestep = push!(Timestep, timestep) From cbf30734b44f97724f190af3d33407fe7242169f Mon Sep 17 00:00:00 2001 From: bkaus Date: Mon, 10 Apr 2023 15:53:32 +0200 Subject: [PATCH 18/18] remove debug messages --- src/read_timestep.jl | 4 ---- test/read_timestep.jl | 1 - 2 files changed, 5 deletions(-) diff --git a/src/read_timestep.jl b/src/read_timestep.jl index ce3bf672..62e18f36 100644 --- a/src/read_timestep.jl +++ b/src/read_timestep.jl @@ -128,7 +128,6 @@ end # The FileName can contain a directory as well; deal with that here function split_path_name(DirName_base::String, FileName::String) FullName = joinpath(DirName_base,FileName) - #id = findlast(Base.Filesystem.path_separator, FullName)[1]; id = findlast("/", FullName)[1]; DirName = FullName[1:id-1] @@ -150,9 +149,7 @@ It will return `data_output` which is a `CartData` output structure. function Read_LaMEM_PVTR_File(DirName_base::String, FileName::String; fields=nothing) CurDir = pwd(); - @show DirName_base FileName DirName, File = split_path_name(DirName_base, FileName) - @show DirName File # change to directory cd(DirName) @@ -237,7 +234,6 @@ function readPVD(FileName::String) # retrieve the timestep file_name = split(file,Base.Filesystem.path_separator)[1]; - #file_name = split(file,"/")[1]; timestep = parse(Int64,split(file_name,"_")[2]); Timestep = push!(Timestep, timestep) diff --git a/test/read_timestep.jl b/test/read_timestep.jl index 9129625b..b36806a5 100644 --- a/test/read_timestep.jl +++ b/test/read_timestep.jl @@ -6,7 +6,6 @@ using LaMEM # Read a timestep FileName="FB_multigrid" Timestep = 1 - @show readdir() data, time = Read_LaMEM_timestep(FileName,Timestep) @test sum(data.fields.phase) ≈ 736.36414f0 @test sum(data.fields.strain_rate[1,:,:,:]) ≈ -0.019376338f0