From cd77f6721f5e60013338720c46046a7155095907 Mon Sep 17 00:00:00 2001 From: Robert Weinbrenner Date: Thu, 11 Dec 2025 10:50:01 +0100 Subject: [PATCH 1/8] Erste Version --- docs/src/howto_parameter-study.md | 115 ++++++++++++++++++++++++++++++ 1 file changed, 115 insertions(+) create mode 100644 docs/src/howto_parameter-study.md diff --git a/docs/src/howto_parameter-study.md b/docs/src/howto_parameter-study.md new file mode 100644 index 00000000..7ec70557 --- /dev/null +++ b/docs/src/howto_parameter-study.md @@ -0,0 +1,115 @@ +# Perform a parameter study + +The following section explains how to perform a parameter study. + +### 1. Define the simulation conditions + +Here is an example function that creates a peridynamic body and defines the boundary conditions of each simulation. +This simulation is intended to simulate wave propagation in the defined body. + +```bash +function job_xwave(setup, root) + (; E, nu, rho, T, vmax, npyz, m) = setup + lx, lyz = 0.2, 0.002 + Δx = lyz / npyz + pos, vol = uniform_box(lx, lyz, lyz, Δx) + ids = sortperm(pos[1,:]) + body = Body(BBMaterial(), pos[:, ids], vol[ids]) + horizon = m * Δx + material!(body; horizon, rho, E, nu) + point_set!(x -> x < -lx / 2 + 1.2Δx, body, :left) + velocity_bc!(t -> t < T ? vmax * sin(2π / T * t) : 0.0, body, :left, :x) + vv = VelocityVerlet(time=1e-4, safety_factor=0.7) + simname = @sprintf("xwave_npyz-%d_m-%d", npyz, round(Int, m)) + path = joinpath(root, simname) + job = Job(body, vv; path, freq=5, fields=(:displacement,)) + return job +end +``` +The line +```bash +(; E, nu, rho, T, vmax, npyz, m) = setup +``` +creates a setup in which the material parameters are defined. + + +Now, various setups are created in which the respective parameters are defined as desired. A separate simulation is performed for each setup. Here the parameter `m`, which controls the peridynamic horizon `h = m * Δx`, varies. + +```bash +setups = [ + (; E=210e9, nu=0.25, rho=7850.0, T=1.0e-5, vmax=2.0, npyz=4, m=3.015), + (; E=210e9, nu=0.25, rho=7850.0, T=1.0e-5, vmax=2.0, npyz=4, m=4.015), + (; E=210e9, nu=0.25, rho=7850.0, T=1.0e-5, vmax=2.0, npyz=4, m=5.015), +] +``` + +The location where the simulation results are to be stored can be specified with `joinpath` +```bash +root = joinpath("results", "xwave_study") +``` + +### 2. Submit study + +To execute the previously defined setups one after the other, first create a study object and then call `submit!`: +``` +study = Study(job_xwave, setups; root) +submit!(study) +``` + + +### 3. Postprocessing + +Now that the simulations are complete, post-processing can begin. This is necessary in order to calculate and display the actual values needed from the raw data generated. In this example the needed values are the theoretical wave velocity `c_0`, the measured wave velocity from the simulation `c_w` and the difference between them two `dc` as well as the percentage deviation `dcp`. + +#### 3.1 Postprocessing function + +First, a standard result structure for post-processing results is created. Then the `process_each_job` function iterates over all simulations (study) and executes the postprocessing code shown below for each individual setup. + +```bash +default_result = (; c_0=NaN, c_w=NaN, dc=NaN, dcp=NaN) +results = process_each_job(study, default_result) do job, setup +``` +#### 3.2 Create an output file and write the results to it + +To create a CSV file (here `wavepos.csv` resp. `wave_position_data_file`) in a desired path where the data can be saved, the following can be executed. +```bash +post_path = joinpath(job.options.root, "post") +mkpath(post_path) +wave_position_data_file = joinpath(post_path, "wavepos.csv") +``` + +The relevant simulation data (here the values for `t`, `x̂` and `û`) is now saved in the CSV file. +```bash +open(wave_position_data_file, "a+") do io + @printf(io, "%.12f,%.12f,%.12f\n", t, x̂, û) +``` + + +#### 3.3 Calculate and display the values + +With the `readdlm` command the created CSV file is read. Then the resulting columns are assigned to `t`, `x̂`, and `û`. + +```bash +results = readdlm(wave_position_data_file, ',', Float64) +t, x̂, û = results[:,1], results[:,2], results[:,3] +``` + + +Now the actual values can be calculated: + +```bash +c_0 = sqrt(E / rho) +c_w = calc_velocity(t, x̂, û) +dc = c_w - c_0 +dcp = 100 * dc / c_0 +``` + + +Finally, the four determined values are displayed in the terminal: + +```bash +printstyled("\nRESULTS FOR SIMULATION: `$(basename(job.options.root))`:\n", bold=true) +printstyled(@sprintf("c_0 = %8.2f m/s\n", c_0); color=:green, bold=true) +printstyled(@sprintf("c_w = %8.2f m/s\n", c_w); color=:blue, bold=true) +printstyled(@sprintf("Δc = %8.2f m/s (%.3f %%)\n", dc, dcp); color=:red, bold=true) +``` From f79dd2b1623a67ec2d200f842512b2c3e713cc4b Mon Sep 17 00:00:00 2001 From: Robert Weinbrenner Date: Thu, 11 Dec 2025 11:08:07 +0100 Subject: [PATCH 2/8] Erste Version neu --- docs/src/Parameter-study.jl | 93 +++++++++++++++++++++++++++++++++++++ 1 file changed, 93 insertions(+) create mode 100644 docs/src/Parameter-study.jl diff --git a/docs/src/Parameter-study.jl b/docs/src/Parameter-study.jl new file mode 100644 index 00000000..684d7e1d --- /dev/null +++ b/docs/src/Parameter-study.jl @@ -0,0 +1,93 @@ +using Peridynamics +using Printf +using DelimitedFiles + +## PARAMETER STUDY SIMULATIONS +function job_xwave(setup, root) + (; E, nu, rho, T, vmax, npyz, m) = setup + lx, lyz = 0.2, 0.002 + Δx = lyz / npyz + pos, vol = uniform_box(lx, lyz, lyz, Δx) + ids = sortperm(pos[1,:]) + body = Body(BBMaterial(), pos[:, ids], vol[ids]) + horizon = m * Δx + material!(body; horizon, rho, E, nu) + point_set!(x -> x < -lx / 2 + 1.2Δx, body, :left) + velocity_bc!(t -> t < T ? vmax * sin(2π / T * t) : 0.0, body, :left, :x) + vv = VelocityVerlet(time=1e-4, safety_factor=0.7) + simname = @sprintf("xwave_npyz-%d_m-%d", npyz, round(Int, m)) + path = joinpath(root, simname) + job = Job(body, vv; path, freq=5, fields=(:displacement,)) + return job +end +setups = [ + (; E=210e9, nu=0.25, rho=7850.0, T=1.0e-5, vmax=2.0, npyz=4, m=3.015), + (; E=210e9, nu=0.25, rho=7850.0, T=1.0e-5, vmax=2.0, npyz=4, m=4.015), + (; E=210e9, nu=0.25, rho=7850.0, T=1.0e-5, vmax=2.0, npyz=4, m=5.015), +] +root = joinpath("results", "xwave_study") +rm(root; recursive=true, force=true) + + +study = Study(job_xwave, setups; root) + +## SUBMIT STUDY +submit!(study) + +## POST-PROCESSING +function calc_velocity(t, x_w, u_w) # Berechnung der Wellengeschwindigkeit anhand der Position und Verschiebung der Welle + u_0 = first(u_w) # Anfangsverschiebung der Welle + valid_until = findfirst(x -> !isapprox(u_0, x; rtol=0.01), u_w) + if !isnothing(valid_until) # Überprüfung, ob ein gültiger Punkt gefunden wurde + x_w = x_w[1:valid_until-1] # Beschränkung der Positions-Daten auf gültige Punkte + t = t[1:valid_until-1] # Beschränkung der Zeit-Daten auf gültige Punkte + end + n = length(t) # Anzahl der Datenpunkte + @assert n == length(x_w) # Sicherstellen, dass die Längen übereinstimmen + t̄ = sum(t) / n #Berechnung des Mittelwerts der Zeit + x̄ = sum(x_w) / n #Berechnung des Mittelwerts der Position + v = sum((t .- t̄) .* (x_w .- x̄)) / sum((t .- t̄) .^ 2) #Berechnung der Steigung der Regressionsgeraden = Geschwindigkeit + return v # Rückgabe der berechneten Geschwindigkeit +end + +default_result = (; c_0=NaN, c_w=NaN, dc=NaN, dcp=NaN) # Standard-Ergebnisstruktur für Post-Processing-Ergebnisse +results = process_each_job(study, default_result) do job, setup # Post-Processing für jede Simulation im Parameter-Studium + # extract parameters that are needed for post-processing + (; E, rho, T, npyz) = setup + Δx = 0.002 / npyz + + # prepare some output files + post_path = joinpath(job.options.root, "post") # Pfad für Post-Processing-Daten + mkpath(post_path) # required to ensure the directory exists (Erstellt das Verzeichnis für die Post-Processing-Daten) + wave_position_data_file = joinpath(post_path, "wavepos.csv") # Datei zur Speicherung der Wellenpositionsdaten + + # get the wave position over time + process_each_export(job; serial=true) do r0, r, id # Liest die Position und Verschiebung der Welle aus den Export-Dateien + t = first(r[:time]) + t < 1.5T && return nothing # nur Zeitpunkte größer als 1.5T berücksichtigen + ymax = @views maximum(r0[:position][2, :]) # maximale y-position + pids = findall(eachcol(r0[:position])) do point + isapprox(point[2], ymax; atol=1.01Δx/2) + end + û, pids_id = @views findmax(r[:displacement][1, pids]) # maximale Verschiebung der Welle in x-Richtung + x̂ = r[:position][1, pids[pids_id]] # Position der Welle in x-Richtung zum gegebenen Zeitpunkt + open(wave_position_data_file, "a+") do io # speichert die Zeit, Position und Verschiebung der Welle in CSV-Datei + @printf(io, "%.12f,%.12f,%.12f\n", t, x̂, û) + end + return nothing + end + + # calculate wave speed + results = readdlm(wave_position_data_file, ',', Float64) + t, x̂, û = results[:,1], results[:,2], results[:,3] + c_0 = sqrt(E / rho) # theoretische Wellengeschwindigkeit + c_w = calc_velocity(t, x̂, û) # Berechnung der Wellengeschwindigkeit anhand Simulationsdaten + dc = c_w - c_0 # Differenz der Wellengeschwindigkeiten + dcp = 100 * dc / c_0 # prozentuale Abweichung der Wellengeschwindigkeiten + printstyled("\nRESULTS FOR SIMULATION: `$(basename(job.options.root))`:\n", bold=true) + printstyled(@sprintf("c_0 = %8.2f m/s\n", c_0); color=:green, bold=true) + printstyled(@sprintf("c_w = %8.2f m/s\n", c_w); color=:blue, bold=true) + printstyled(@sprintf("Δc = %8.2f m/s (%.3f %%)\n", dc, dcp); color=:red, bold=true) + + return (; c_0, c_w, dc, dcp) +end From eab02ffe64529581709462345f5354bb4e8550a0 Mon Sep 17 00:00:00 2001 From: Robert Weinbrenner Date: Thu, 11 Dec 2025 11:11:54 +0100 Subject: [PATCH 3/8] Erste Version neu --- docs/make.jl | 1 + docs/src/howto_parameter-study.md | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index c3312c0d..4146f890 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -72,6 +72,7 @@ makedocs(; "How-to guides" => [ "howto_mpi.md", "howto_visualization.md", + "howto_parameter-study.md", ], "Tutorials" => [ joinpath("generated", "tutorial_tension_static.md"), diff --git a/docs/src/howto_parameter-study.md b/docs/src/howto_parameter-study.md index 7ec70557..dcddb97c 100644 --- a/docs/src/howto_parameter-study.md +++ b/docs/src/howto_parameter-study.md @@ -33,7 +33,7 @@ The line creates a setup in which the material parameters are defined. -Now, various setups are created in which the respective parameters are defined as desired. A separate simulation is performed for each setup. Here the parameter `m`, which controls the peridynamic horizon `h = m * Δx`, varies. +Now various setups are created in which the respective parameters are defined as desired. A separate simulation is performed for each setup. Here the parameter `m`, which controls the peridynamic horizon `h = m * Δx`, varies. ```bash setups = [ From a19186a5c452beef8eee1b48758a3412f5518908 Mon Sep 17 00:00:00 2001 From: robertweinbrenner Date: Thu, 11 Dec 2025 12:25:29 +0100 Subject: [PATCH 4/8] Update docs/src/howto_parameter-study.md Co-authored-by: Kai Partmann <68582683+kaipartmann@users.noreply.github.com> --- docs/src/howto_parameter-study.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/howto_parameter-study.md b/docs/src/howto_parameter-study.md index dcddb97c..23c12a70 100644 --- a/docs/src/howto_parameter-study.md +++ b/docs/src/howto_parameter-study.md @@ -89,7 +89,7 @@ open(wave_position_data_file, "a+") do io With the `readdlm` command the created CSV file is read. Then the resulting columns are assigned to `t`, `x̂`, and `û`. -```bash +```julia results = readdlm(wave_position_data_file, ',', Float64) t, x̂, û = results[:,1], results[:,2], results[:,3] ``` From 0db131f7927bd35a79051c8fa00113fd650e9dde Mon Sep 17 00:00:00 2001 From: robertweinbrenner Date: Thu, 11 Dec 2025 13:05:39 +0100 Subject: [PATCH 5/8] Update docs/src/howto_parameter-study.md Co-authored-by: Kai Partmann <68582683+kaipartmann@users.noreply.github.com> --- docs/src/howto_parameter-study.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/src/howto_parameter-study.md b/docs/src/howto_parameter-study.md index 23c12a70..984471b2 100644 --- a/docs/src/howto_parameter-study.md +++ b/docs/src/howto_parameter-study.md @@ -79,9 +79,10 @@ wave_position_data_file = joinpath(post_path, "wavepos.csv") ``` The relevant simulation data (here the values for `t`, `x̂` and `û`) is now saved in the CSV file. -```bash +```julia open(wave_position_data_file, "a+") do io @printf(io, "%.12f,%.12f,%.12f\n", t, x̂, û) +end ``` From c1234dca7be4fcc2fe67aaae54d1394f631a7fc1 Mon Sep 17 00:00:00 2001 From: Robert Weinbrenner Date: Wed, 14 Jan 2026 10:15:24 +0100 Subject: [PATCH 6/8] Update regarding the comments --- docs/src/howto_parameter-study.md | 72 +++++++++++++++++++++++++------ 1 file changed, 59 insertions(+), 13 deletions(-) diff --git a/docs/src/howto_parameter-study.md b/docs/src/howto_parameter-study.md index 984471b2..3a0454ee 100644 --- a/docs/src/howto_parameter-study.md +++ b/docs/src/howto_parameter-study.md @@ -7,7 +7,7 @@ The following section explains how to perform a parameter study. Here is an example function that creates a peridynamic body and defines the boundary conditions of each simulation. This simulation is intended to simulate wave propagation in the defined body. -```bash +```julia function job_xwave(setup, root) (; E, nu, rho, T, vmax, npyz, m) = setup lx, lyz = 0.2, 0.002 @@ -26,8 +26,9 @@ function job_xwave(setup, root) return job end ``` + The line -```bash +```julia (; E, nu, rho, T, vmax, npyz, m) = setup ``` creates a setup in which the material parameters are defined. @@ -35,7 +36,7 @@ creates a setup in which the material parameters are defined. Now various setups are created in which the respective parameters are defined as desired. A separate simulation is performed for each setup. Here the parameter `m`, which controls the peridynamic horizon `h = m * Δx`, varies. -```bash +```julia setups = [ (; E=210e9, nu=0.25, rho=7850.0, T=1.0e-5, vmax=2.0, npyz=4, m=3.015), (; E=210e9, nu=0.25, rho=7850.0, T=1.0e-5, vmax=2.0, npyz=4, m=4.015), @@ -44,14 +45,14 @@ setups = [ ``` The location where the simulation results are to be stored can be specified with `joinpath` -```bash +```julia root = joinpath("results", "xwave_study") ``` ### 2. Submit study To execute the previously defined setups one after the other, first create a study object and then call `submit!`: -``` +```julia study = Study(job_xwave, setups; root) submit!(study) ``` @@ -59,29 +60,76 @@ submit!(study) ### 3. Postprocessing -Now that the simulations are complete, post-processing can begin. This is necessary in order to calculate and display the actual values needed from the raw data generated. In this example the needed values are the theoretical wave velocity `c_0`, the measured wave velocity from the simulation `c_w` and the difference between them two `dc` as well as the percentage deviation `dcp`. +Now that the simulations are complete, post-processing can begin. This is necessary in order to calculate and display the actual values needed from the raw data generated. In this example the needed values are the theoretical wave velocity `c_0`, the measured wave velocity from the simulation `c_w` and the difference between them two `dc` as well as the percentage deviation `dcp`. + +The complete postprocessing code is shown below and will be discussed in the following sections. + + + +```julia +default_result = (; c_0=NaN, c_w=NaN, dc=NaN, dcp=NaN) +results = process_each_job(study, default_result) do job, setup + (; E, rho, T, npyz) = setup + Δx = 0.002 / npyz + + post_path = joinpath(job.options.root, "post") + mkpath(post_path) + wave_position_data_file = joinpath(post_path, "wavepos.csv") + + process_each_export(job; serial=true) do r0, r, id + t = first(r[:time]) + t < 1.5T && return nothing + ymax = @views maximum(r0[:position][2, :]) + pids = findall(eachcol(r0[:position])) do point + isapprox(point[2], ymax; atol=1.01Δx/2) + end + û, pids_id = @views findmax(r[:displacement][1, pids]) + x̂ = r[:position][1, pids[pids_id]] + open(wave_position_data_file, "a+") do io + @printf(io, "%.12f,%.12f,%.12f\n", t, x̂, û) + end + return nothing + end + + results = readdlm(wave_position_data_file, ',', Float64) + t, x̂, û = results[:,1], results[:,2], results[:,3] + c_0 = sqrt(E / rho) + c_w = calc_velocity(t, x̂, û) + dc = c_w - c_0 + dcp = 100 * dc / c_0 + printstyled("\nRESULTS FOR SIMULATION: `$(basename(job.options.root))`:\n", bold=true) + printstyled(@sprintf("c_0 = %8.2f m/s\n", c_0); color=:green, bold=true) + printstyled(@sprintf("c_w = %8.2f m/s\n", c_w); color=:blue, bold=true) + printstyled(@sprintf("Δc = %8.2f m/s (%.3f %%)\n", dc, dcp); color=:red, bold=true) + + return (; c_0, c_w, dc, dcp) +end +``` #### 3.1 Postprocessing function First, a standard result structure for post-processing results is created. Then the `process_each_job` function iterates over all simulations (study) and executes the postprocessing code shown below for each individual setup. -```bash +```julia default_result = (; c_0=NaN, c_w=NaN, dc=NaN, dcp=NaN) results = process_each_job(study, default_result) do job, setup ``` + + #### 3.2 Create an output file and write the results to it To create a CSV file (here `wavepos.csv` resp. `wave_position_data_file`) in a desired path where the data can be saved, the following can be executed. -```bash +```julia post_path = joinpath(job.options.root, "post") mkpath(post_path) wave_position_data_file = joinpath(post_path, "wavepos.csv") ``` -The relevant simulation data (here the values for `t`, `x̂` and `û`) is now saved in the CSV file. ```julia open(wave_position_data_file, "a+") do io @printf(io, "%.12f,%.12f,%.12f\n", t, x̂, û) +end + @printf(io, "%.12f,%.12f,%.12f\n", t, x̂, û) end ``` @@ -97,8 +145,7 @@ t, x̂, û = results[:,1], results[:,2], results[:,3] Now the actual values can be calculated: - -```bash +```julia c_0 = sqrt(E / rho) c_w = calc_velocity(t, x̂, û) dc = c_w - c_0 @@ -107,8 +154,7 @@ dcp = 100 * dc / c_0 Finally, the four determined values are displayed in the terminal: - -```bash +```julia printstyled("\nRESULTS FOR SIMULATION: `$(basename(job.options.root))`:\n", bold=true) printstyled(@sprintf("c_0 = %8.2f m/s\n", c_0); color=:green, bold=true) printstyled(@sprintf("c_w = %8.2f m/s\n", c_w); color=:blue, bold=true) From 1b4f2ec36240c547236cb6682d52daf1c0595a54 Mon Sep 17 00:00:00 2001 From: Kai Partmann Date: Wed, 14 Jan 2026 10:35:14 +0100 Subject: [PATCH 7/8] Minor hints regarding the do-syntax --- docs/src/howto_parameter-study.md | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/docs/src/howto_parameter-study.md b/docs/src/howto_parameter-study.md index 3a0454ee..b738281e 100644 --- a/docs/src/howto_parameter-study.md +++ b/docs/src/howto_parameter-study.md @@ -113,8 +113,17 @@ First, a standard result structure for post-processing results is created. Then ```julia default_result = (; c_0=NaN, c_w=NaN, dc=NaN, dcp=NaN) results = process_each_job(study, default_result) do job, setup + (; E, rho, T, npyz) = setup + Δx = 0.002 / npyz + + # ... postprocessing code for each job (see sections 3.2-3.3) ... + + return (; c_0, c_w, dc, dcp) +end ``` +The `do` block receives `job` and `setup` as parameters and should return the computed results for each simulation. + #### 3.2 Create an output file and write the results to it @@ -128,8 +137,6 @@ wave_position_data_file = joinpath(post_path, "wavepos.csv") ```julia open(wave_position_data_file, "a+") do io @printf(io, "%.12f,%.12f,%.12f\n", t, x̂, û) -end - @printf(io, "%.12f,%.12f,%.12f\n", t, x̂, û) end ``` From e16a2aee771959a637b4e75d6aaf6ae70d874c16 Mon Sep 17 00:00:00 2001 From: Robert Weinbrenner Date: Thu, 12 Feb 2026 11:09:43 +0100 Subject: [PATCH 8/8] Update of the postprocessing part --- docs/src/howto_parameter-study.md | 119 ++++++++++++++++-------------- 1 file changed, 65 insertions(+), 54 deletions(-) diff --git a/docs/src/howto_parameter-study.md b/docs/src/howto_parameter-study.md index b738281e..c8d36f41 100644 --- a/docs/src/howto_parameter-study.md +++ b/docs/src/howto_parameter-study.md @@ -60,55 +60,58 @@ submit!(study) ### 3. Postprocessing -Now that the simulations are complete, post-processing can begin. This is necessary in order to calculate and display the actual values needed from the raw data generated. In this example the needed values are the theoretical wave velocity `c_0`, the measured wave velocity from the simulation `c_w` and the difference between them two `dc` as well as the percentage deviation `dcp`. +Now that the simulations are complete, the post-processing can begin. This is necessary in order to calculate and display the actual values needed from the raw data generated. In this example the needed values are the theoretical wave velocity `c_0`, the measured wave velocity from the simulation `c_w` and the difference between them two `dc` as well as the percentage deviation `dcp`. The complete postprocessing code is shown below and will be discussed in the following sections. - - ```julia default_result = (; c_0=NaN, c_w=NaN, dc=NaN, dcp=NaN) results = process_each_job(study, default_result) do job, setup (; E, rho, T, npyz) = setup Δx = 0.002 / npyz - post_path = joinpath(job.options.root, "post") - mkpath(post_path) - wave_position_data_file = joinpath(post_path, "wavepos.csv") - - process_each_export(job; serial=true) do r0, r, id + default_exp_ret = (; time=NaN, wave_position=NaN, amplitude_ux=NaN) + res = process_each_export(job, default_exp_ret) do r0, r, id t = first(r[:time]) - t < 1.5T && return nothing + t < 1.5T && return default_exp_ret ymax = @views maximum(r0[:position][2, :]) pids = findall(eachcol(r0[:position])) do point isapprox(point[2], ymax; atol=1.01Δx/2) end û, pids_id = @views findmax(r[:displacement][1, pids]) x̂ = r[:position][1, pids[pids_id]] - open(wave_position_data_file, "a+") do io - @printf(io, "%.12f,%.12f,%.12f\n", t, x̂, û) - end - return nothing + return (; time=t, wave_position=x̂, amplitude_ux=û) end - results = readdlm(wave_position_data_file, ',', Float64) - t, x̂, û = results[:,1], results[:,2], results[:,3] - c_0 = sqrt(E / rho) - c_w = calc_velocity(t, x̂, û) - dc = c_w - c_0 - dcp = 100 * dc / c_0 - printstyled("\nRESULTS FOR SIMULATION: `$(basename(job.options.root))`:\n", bold=true) - printstyled(@sprintf("c_0 = %8.2f m/s\n", c_0); color=:green, bold=true) - printstyled(@sprintf("c_w = %8.2f m/s\n", c_w); color=:blue, bold=true) - printstyled(@sprintf("Δc = %8.2f m/s (%.3f %%)\n", dc, dcp); color=:red, bold=true) - - return (; c_0, c_w, dc, dcp) -end + allres = (; (k => getfield.(res, k) for k in keys(default_exp_ret))...) + (; time, wave_position, amplitude_ux) = allres + idxs = findall(x -> !isnan(x), time) + + c_0 = sqrt(E / rho) + c_w = calc_velocity(time[idxs], wave_position[idxs], amplitude_ux[idxs]) + dc = c_w - c_0 + dcp = 100 * dc / c_0 + + @mpiroot begin + bold = true + simname = basename(job.options.root) + printstyled("\nRESULTS FOR SIMULATION: `$(simname)`:\n"; bold) + printstyled(@sprintf("c_0 = %8.2f m/s\n", c_0); color=:green, bold) + printstyled(@sprintf("c_w = %8.2f m/s\n", c_w); color=:blue, bold) + printstyled(@sprintf("Δc = %8.2f m/s (%.3f %%)\n", dc, dcp); color=:red, bold) + end + + return (; c_0, c_w, dc, dcp) +end ``` + #### 3.1 Postprocessing function -First, a standard result structure for post-processing results is created. Then the `process_each_job` function iterates over all simulations (study) and executes the postprocessing code shown below for each individual setup. +First, `default_result` is used to define a default result structure (a NamedTuple). +The function `process_each_job` then runs through all simulations (`study`) and executes the postprocessing code listed in sections 3.2-3 for each individual setup. The results are stored in `res`. +The `do` block receives `job` and `setup` as parameters and returns the calculated results for each simulation. + ```julia default_result = (; c_0=NaN, c_w=NaN, dc=NaN, dcp=NaN) @@ -122,48 +125,56 @@ results = process_each_job(study, default_result) do job, setup end ``` -The `do` block receives `job` and `setup` as parameters and should return the computed results for each simulation. +#### 3.2 Identify irrelevant data -#### 3.2 Create an output file and write the results to it - -To create a CSV file (here `wavepos.csv` resp. `wave_position_data_file`) in a desired path where the data can be saved, the following can be executed. -```julia -post_path = joinpath(job.options.root, "post") -mkpath(post_path) -wave_position_data_file = joinpath(post_path, "wavepos.csv") -``` +Now, invalid or irrelevant data that is not to be processed further is collected. +For this purpose, another named tuple `default_exp_ret` is defined as the result structure. The function `process_each_export` iterates over the stored export data of the job and then stores only the relevant results in `res`. ```julia -open(wave_position_data_file, "a+") do io - @printf(io, "%.12f,%.12f,%.12f\n", t, x̂, û) +default_exp_ret = (; time=NaN, wave_position=NaN, amplitude_ux=NaN) +res = process_each_export(job, default_exp_ret) do r0, r, id + t = first(r[:time]) + t < 1.5T && return default_exp_ret + ymax = @views maximum(r0[:position][2, :]) + pids = findall(eachcol(r0[:position])) do point + isapprox(point[2], ymax; atol=1.01Δx/2) + end + û, pids_id = @views findmax(r[:displacement][1, pids]) + x̂ = r[:position][1, pids[pids_id]] + return (; time=t, wave_position=x̂, amplitude_ux=û) end ``` -#### 3.3 Calculate and display the values +#### 3.3 Transformation of the results -With the `readdlm` command the created CSV file is read. Then the resulting columns are assigned to `t`, `x̂`, and `û`. +Now, from the results stored in `res`, the vectors of the fields (`time`, `wave_position`, `amplitude_ux`) are collected in a NamedTuple, unpacked into local vectors, and the indices of the non-NaN entries are determined with `findall` so that only valid measurement points are used later. ```julia -results = readdlm(wave_position_data_file, ',', Float64) -t, x̂, û = results[:,1], results[:,2], results[:,3] +allres = (; (k => getfield.(res, k) for k in keys(default_exp_ret))...) +(; time, wave_position, amplitude_ux) = allres +idxs = findall(x -> !isnan(x), time) ``` -Now the actual values can be calculated: -```julia -c_0 = sqrt(E / rho) -c_w = calc_velocity(t, x̂, û) -dc = c_w - c_0 -dcp = 100 * dc / c_0 -``` +#### 3.4 Calculate and display the values + +Now the actual values can be calculated and displayed: -Finally, the four determined values are displayed in the terminal: ```julia -printstyled("\nRESULTS FOR SIMULATION: `$(basename(job.options.root))`:\n", bold=true) -printstyled(@sprintf("c_0 = %8.2f m/s\n", c_0); color=:green, bold=true) -printstyled(@sprintf("c_w = %8.2f m/s\n", c_w); color=:blue, bold=true) -printstyled(@sprintf("Δc = %8.2f m/s (%.3f %%)\n", dc, dcp); color=:red, bold=true) + c_0 = sqrt(E / rho) + c_w = calc_velocity(time[idxs], wave_position[idxs], amplitude_ux[idxs]) + dc = c_w - c_0 + dcp = 100 * dc / c_0 + + @mpiroot begin + bold = true + simname = basename(job.options.root) + printstyled("\nRESULTS FOR SIMULATION: `$(simname)`:\n"; bold) + printstyled(@sprintf("c_0 = %8.2f m/s\n", c_0); color=:green, bold) + printstyled(@sprintf("c_w = %8.2f m/s\n", c_w); color=:blue, bold) + printstyled(@sprintf("Δc = %8.2f m/s (%.3f %%)\n", dc, dcp); color=:red, bold) + end ```