From db21a64b3a8dd13a8b6fecc75e28a0e3b49f9b29 Mon Sep 17 00:00:00 2001 From: David Eldon Date: Wed, 4 Oct 2023 14:04:37 -0700 Subject: [PATCH 01/12] Make geqdsk_to_imas modify in place --- src/SD4SOLPS.jl | 10 +++++----- test/runtests.jl | 10 +++++----- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/SD4SOLPS.jl b/src/SD4SOLPS.jl index 2fbeb29..11b032a 100644 --- a/src/SD4SOLPS.jl +++ b/src/SD4SOLPS.jl @@ -7,7 +7,7 @@ using Interpolations: Interpolations #import GGDUtils export find_files_in_allowed_folders -export geqdsk_to_imas +export geqdsk_to_imas! include("$(@__DIR__)/supersize_profile.jl") include("$(@__DIR__)/repair_eq.jl") @@ -63,12 +63,12 @@ function find_files_in_allowed_folders( end """ - geqdsk_to_imas() + geqdsk_to_imas!() Transfers the equilibrium reconstruction in an EFIT-style gEQDSK file into the IMAS DD structure. """ -function geqdsk_to_imas(eqdsk_file, dd; time_index=1) +function geqdsk_to_imas!(eqdsk_file, dd; time_index=1) # https://github.com/JuliaFusion/EFIT.jl/blob/master/src/io.jl g = EFIT.readg(eqdsk_file) # Copying ideas from OMFIT: omfit/omfit_classes/omfit_eqdsk.py / to_omas() @@ -134,7 +134,7 @@ function geqdsk_to_imas(eqdsk_file, dd; time_index=1) limiter.type.description = "first wall" resize!(limiter.unit, 1) limiter.unit[1].outline.r = g.rlim - return limiter.unit[1].outline.z = g.zlim + limiter.unit[1].outline.z = g.zlim end """ @@ -238,7 +238,7 @@ function preparation( println(" eqdsk = ", eqdsk) dd = SOLPS2IMAS.solps2imas(b2fgmtry, b2time, gridspec, b2mn) - geqdsk_to_imas(eqdsk, dd) + geqdsk_to_imas!(eqdsk, dd) println("Loaded input data into IMAS DD") fill_in_extrapolated_core_profile!(dd, "electrons.density"; method=core_method) diff --git a/test/runtests.jl b/test/runtests.jl index a8eebcf..d10577b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -185,7 +185,7 @@ if args["core_profile_extension"] @test isfile(gridspec) @test isfile(eqdsk) dd = SOLPS2IMAS.solps2imas(b2fgmtry, b2time, gridspec, b2mn) - SD4SOLPS.geqdsk_to_imas(eqdsk, dd) + SD4SOLPS.geqdsk_to_imas!(eqdsk, dd) rho = dd.equilibrium.time_slice[1].profiles_1d.rho_tor_norm if !SD4SOLPS.check_rho_1d(dd; time_slice=1) @@ -221,7 +221,7 @@ if args["edge_profile_extension"] # Test for getting mesh spacing b2fgmtry, b2time, b2mn, gridspec, eqdsk = define_default_sample_set() dd = SOLPS2IMAS.solps2imas(b2fgmtry, b2time, gridspec, b2mn) - SD4SOLPS.geqdsk_to_imas(eqdsk, dd) + SD4SOLPS.geqdsk_to_imas!(eqdsk, dd) dpsin = SD4SOLPS.mesh_psi_spacing(dd) @test dpsin > 0.0 @@ -260,7 +260,7 @@ if args["heavy_utilities"] dd = OMAS.dd() eqdsk_file = splitdir(pathof(SD4SOLPS))[1] * "/../sample/geqdsk_iter_small_sample" - SD4SOLPS.geqdsk_to_imas(eqdsk_file, dd) + SD4SOLPS.geqdsk_to_imas!(eqdsk_file, dd) quantity = "electrons.density" prof_time_idx = eq_time_idx = 1 resize!(dd.core_profiles.profiles_1d, prof_time_idx) @@ -293,7 +293,7 @@ if args["repair_eq"] # Prepare sample dd = OMAS.dd() eqdsk = splitdir(pathof(SD4SOLPS))[1] * "/../sample/geqdsk_iter_small_sample" - SD4SOLPS.geqdsk_to_imas(eqdsk, dd) + SD4SOLPS.geqdsk_to_imas!(eqdsk, dd) # Make sure rho is missing nt = length(dd.equilibrium.time_slice) for it ∈ 1:nt @@ -324,7 +324,7 @@ if args["geqdsk_to_imas"] for sample_file ∈ sample_files println(sample_file) dd = OMAS.dd() - SD4SOLPS.geqdsk_to_imas(sample_file, dd; time_index=tslice) + SD4SOLPS.geqdsk_to_imas!(sample_file, dd; time_index=tslice) eqt = dd.equilibrium.time_slice[tslice] # global From 46d6a57aaf4eec5397fb11fd82eda62973ebb94f Mon Sep 17 00:00:00 2001 From: David Eldon Date: Thu, 2 Nov 2023 11:55:53 -0700 Subject: [PATCH 02/12] bugfix in repair_eq --- src/repair_eq.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/repair_eq.jl b/src/repair_eq.jl index d90b293..24b53bb 100644 --- a/src/repair_eq.jl +++ b/src/repair_eq.jl @@ -48,7 +48,7 @@ function add_rho_to_equilibrium!(dd::OMAS.dd) psi = eqt.profiles_1d.psi n = length(psi) if (length(eqt.profiles_1d.rho_tor_norm) > 0) - if max(eqt.profiles_1d.rho_tor_norm) > 0 + if maximum(eqt.profiles_1d.rho_tor_norm) > 0 println("Slice #", it, " already has a rho_tor_norm profile; skipping") continue end From 5a06a39867e4bcc31329e3989b6af9cdc93f43a0 Mon Sep 17 00:00:00 2001 From: David Eldon Date: Wed, 1 Nov 2023 14:54:08 -0700 Subject: [PATCH 03/12] geqdsk to IMAS function records X-points to DD after EFIT finds them --- src/SD4SOLPS.jl | 33 +++++++++++++++++++++++++++++++++ test/runtests.jl | 10 ++++++++++ 2 files changed, 43 insertions(+) diff --git a/src/SD4SOLPS.jl b/src/SD4SOLPS.jl index 11b032a..84a55f4 100644 --- a/src/SD4SOLPS.jl +++ b/src/SD4SOLPS.jl @@ -122,6 +122,39 @@ function geqdsk_to_imas!(eqdsk_file, dd; time_index=1) gq.q_min.rho_tor_norm = g.rhovn[qmin_idx] end + # X-points + xrs, xzs, xpsins, xseps = EFIT.x_points(g) + if length(xrs) > 0 + bx = eqt.boundary.x_point + resize!(bx, length(xrs)) + for i ∈ length(xrs) + bx[i].r = xrs[i] + bx[i].z = xzs[i] + end + nprim = sum(xseps .== 1) + if nprim > 0 + bsx = eqt.boundary_separatrix.x_point + resize!(bsx, nprim) + xrprim = xrs[xseps .== 1] + xzprim = xzs[xseps .== 1] + for i ∈ nprim + bsx[i].r = xrprim[i] + bsx[i].z = xzprim[i] + end + end + nsec = sum(xseps .== 2) + if nsec > 0 + bssx = eqt.boundary_secondary_separatrix.x_point + resize!(bssx, nsec) + xrsec = xrs[xseps .== 2] + xzsec = xzs[xseps .== 2] + for i ∈ nsec + bssx[i].r = xrsec[i] + bssx[i].z = xzsec[i] + end + end + end + # Boundary / LCFS eqt.boundary.outline.r = g.rbbbs eqt.boundary.outline.z = g.zbbbs diff --git a/test/runtests.jl b/test/runtests.jl index d10577b..116fcf8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -376,6 +376,16 @@ if args["geqdsk_to_imas"] # derived @test gq.q_axis == p1.q[1] + # X-points + bx = eqt.boundary.x_point + bsx = eqt.boundary_separatrix.x_point + bssx = eqt.boundary_secondary_separatrix.x_point + nxpt = length(bx) + nprim = length(bsx) + nsec = length(bssx) + @test nxpt >= 1 + @test nxpt >= (nprim + nsec) + # wall limiter = dd.wall.description_2d[1].limiter @test length(limiter.unit[1].outline.r) > 10 From df0d567f80c9f333ce28256b9566077de7bf0527 Mon Sep 17 00:00:00 2001 From: David Eldon Date: Wed, 4 Oct 2023 16:41:46 -0700 Subject: [PATCH 04/12] Add function for creating and recording extended mesh - Find gridlines by starting at inner boundary of existing mesh and following gradients of psi_N. Stop at regular intervals to get nodes in the new mesh. - Modify steepest descent path near secondary X-point - The path that is probably meant to go through the secondary X-point will miss the secondary X-point if there's any numerical error (so, always) - Instead, find the X-point (a proposed change to EFIT.jl does this, but we can move that function back to SD4SOLPS if EFIT.jl doesn't accept it) and draw a line from the secondary X-point to the closest vertex on the inner edge of the mesh, replacing the whacky curved path that almost hits the X-point. - At this point, some of the cells are outside the limiting surface --- Project.toml | 1 + apply_format.sh | 2 + src/SD4SOLPS.jl | 11 +- src/supersize_profile.jl | 390 ++++++++++++++++++++++++++++++++++++++- test/runtests.jl | 39 +++- 5 files changed, 430 insertions(+), 13 deletions(-) create mode 100755 apply_format.sh diff --git a/Project.toml b/Project.toml index 8b2bc92..e1d5140 100644 --- a/Project.toml +++ b/Project.toml @@ -13,6 +13,7 @@ OMAS = "91cfaa06-6526-4804-8666-b540b3feef2f" PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab" PlotUtils = "995b91a9-d308-5afd-9ec6-746e21dbc043" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +PolygonOps = "647866c9-e3ac-4575-94e7-e3d426903924" SOLPS2IMAS = "09becab6-0636-4c23-a92a-2b3723265c31" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" diff --git a/apply_format.sh b/apply_format.sh new file mode 100755 index 0000000..0101490 --- /dev/null +++ b/apply_format.sh @@ -0,0 +1,2 @@ +#!/bin/bash +julia -e 'using JuliaFormatter; format(".", verbose=true)' diff --git a/src/SD4SOLPS.jl b/src/SD4SOLPS.jl index 84a55f4..09a9f8e 100644 --- a/src/SD4SOLPS.jl +++ b/src/SD4SOLPS.jl @@ -123,7 +123,7 @@ function geqdsk_to_imas!(eqdsk_file, dd; time_index=1) end # X-points - xrs, xzs, xpsins, xseps = EFIT.x_points(g) + xrs, xzs, xpsins, xseps = EFIT.x_points(g; within_limiter_only=false) if length(xrs) > 0 bx = eqt.boundary.x_point resize!(bx, length(xrs)) @@ -135,8 +135,8 @@ function geqdsk_to_imas!(eqdsk_file, dd; time_index=1) if nprim > 0 bsx = eqt.boundary_separatrix.x_point resize!(bsx, nprim) - xrprim = xrs[xseps .== 1] - xzprim = xzs[xseps .== 1] + xrprim = xrs[xseps.==1] + xzprim = xzs[xseps.==1] for i ∈ nprim bsx[i].r = xrprim[i] bsx[i].z = xzprim[i] @@ -146,8 +146,8 @@ function geqdsk_to_imas!(eqdsk_file, dd; time_index=1) if nsec > 0 bssx = eqt.boundary_secondary_separatrix.x_point resize!(bssx, nsec) - xrsec = xrs[xseps .== 2] - xzsec = xzs[xseps .== 2] + xrsec = xrs[xseps.==2] + xzsec = xzs[xseps.==2] for i ∈ nsec bssx[i].r = xrsec[i] bssx[i].z = xzsec[i] @@ -168,6 +168,7 @@ function geqdsk_to_imas!(eqdsk_file, dd; time_index=1) resize!(limiter.unit, 1) limiter.unit[1].outline.r = g.rlim limiter.unit[1].outline.z = g.zlim + return end """ diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index 5e8cedc..333f4bc 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -6,10 +6,13 @@ Utilities for extrapolating profiles using OMAS: OMAS using Interpolations: Interpolations using GGDUtils: GGDUtils +using PolygonOps: PolygonOps +using SOLPS2IMAS: SOLPS2IMAS export extrapolate_core export fill_in_extrapolated_core_profile export mesh_psi_spacing +export find_x_points! """ cumul_integrate(x::AbstractVector, y::AbstractVector) @@ -205,10 +208,8 @@ function fill_in_extrapolated_core_profile!( in_bounds = (r .< maximum(r_eq)) .& (r .> minimum(r_eq)) .& (z .> minimum(z_eq)) .& (z .< maximum(z_eq)) - # println(in_bounds) psi_for_quantity = 10.0 .+ zeros(length(r)) psi_for_quantity[in_bounds] = rzpi.(r[in_bounds], z[in_bounds]) - # println(length(psi1_eq), ", ", length(rho1_eq)) rho_for_quantity = copy(psi_for_quantity) in_bounds = psi_for_quantity .<= 1.0 dpsi = diff(psi1_eq) @@ -305,6 +306,10 @@ grid_ggd_idx: index of the grid_ggd to use. For a typical SOLPS run, the SOLPS g is used, then this index will need to be specified. space_idx: index of the space to use. For a typical SOLPS run, there will be only one space so this index will mostly remain at 1. +avoid_guard_cell: assume that the last cell is a guard cell so take end-2 and end-1 + instead of end and end-1 +spacing_rule: "edge" or "mean" to make spacing of new cells (in psi_N) be the same + as the spacing at the edge of the mesh, or the same as the average spacing """ #! format on function mesh_psi_spacing( @@ -313,6 +318,8 @@ function mesh_psi_spacing( eq_profiles_2d_idx::Int64=1, grid_ggd_idx::Int64=1, space_idx::Int64=1, + avoid_guard_cell::Bool=true, + spacing_rule="mean", ) # Inspect input if length(dd.equilibrium.time_slice) < eq_time_idx @@ -344,8 +351,8 @@ function mesh_psi_spacing( psib = eqt.global_quantities.psi_boundary psin_eq = (psi .- psia) ./ (psib - psia) rzpi = Interpolations.linear_interpolation((r_eq, z_eq), psin_eq) - println(minimum(r_eq), ", ", maximum(r_eq)) - println(minimum(z_eq), ", ", maximum(z_eq)) + # println(minimum(r_eq), ", ", maximum(r_eq)) + # println(minimum(z_eq), ", ", maximum(z_eq)) # Get a row of cells. Since the mesh should be aligned to the flux surfaces, # it shouldn't matter which row is used, although the divertor rows might be @@ -357,16 +364,385 @@ function mesh_psi_spacing( midplane_cell_centers = GGDUtils.get_subset_centers(space, midplane_subset) r_mesh = [midplane_cell_centers[i][1] for i ∈ eachindex(midplane_cell_centers)] z_mesh = [midplane_cell_centers[i][2] for i ∈ eachindex(midplane_cell_centers)] - println(minimum(r_mesh), ", ", maximum(r_mesh)) - println(minimum(z_mesh), ", ", maximum(z_mesh)) psin_mesh = rzpi.(r_mesh, z_mesh) # This should come out sorted, but with GGD, who knows. ii = sortperm(psin_mesh) psin = psin_mesh[ii] - dpsin = psin[end] - psin[end-1] + if spacing_rule == "edge" + if avoid_guard_cell + dpsin = psin[end-1] - psin[end-2] + else + dpsin = psin[end] - psin[end-1] + end + else + if avoid_guard_cell + dpsin = diff(psin[2:end-1]) + else + dpsin = diff(psin) + end + dpsin = sum(dpsin) / length(dpsin) + end return dpsin end +""" + function mesh_extension_sol() + +Extends the mesh out into the SOL +""" +function mesh_extension_sol!( + dd::OMAS.dd; + eq_time_idx::Int64=1, + eq_profiles_2d_idx::Int64=1, + grid_ggd_idx::Int64=1, + space_idx::Int64=1, +) + # Prep flux map + eqt = dd.equilibrium.time_slice[eq_time_idx] + p2 = eqt.profiles_2d[eq_profiles_2d_idx] + r_eq = p2.grid.dim1 + z_eq = p2.grid.dim2 + psi = p2.psi + psia = eqt.global_quantities.psi_axis + psib = eqt.global_quantities.psi_boundary + psin_eq = (psi .- psia) ./ (psib - psia) + rzpi = Interpolations.LinearInterpolation((r_eq, z_eq), psin_eq) + + # Use wall to find maximum extent of contouring project + limiter = dd.wall.description_2d[1].limiter + wall_r = limiter.unit[1].outline.r + wall_z = limiter.unit[1].outline.z + wall_psin = rzpi.(wall_r, wall_z) + + # Use ggd mesh to find inner limit of contouring project + grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] + space = grid_ggd.space[space_idx] + midplane_subset = + SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 11) + midplane_cell_centers = GGDUtils.get_subset_centers(space, midplane_subset) + psin_midplane = rzpi.(midplane_cell_centers[end][1], midplane_cell_centers[end][2]) + + # Choose contour levels + # The psi spacing function doesn't do much that's unique anymore. + # Should it be absorbed into this function intead? + dpsin = mesh_psi_spacing( + dd; + eq_time_idx=eq_time_idx, + eq_profiles_2d_idx=eq_profiles_2d_idx, + grid_ggd_idx=grid_ggd_idx, + space_idx=space_idx, + avoid_guard_cell=true, + spacing_rule="mean", + ) + lvlstart = maximum(psin_midplane * sign(dpsin)) / sign(dpsin) + dpsin + lvlend = maximum(wall_psin * sign(dpsin)) / sign(dpsin) + dpsin + nlvl = Int64(ceil((lvlend - lvlstart) / dpsin)) + psin_levels = collect(LinRange(lvlstart, lvlend, nlvl)) + if hasproperty(eqt.boundary_secondary_separatrix, :psi) + secondary_psi = eqt.boundary_secondary_separatrix.psi + secondary_psin = (secondary_psi - psia) / (psib - psia) + else + secondary_psin = lvlend + dpsin + end + + # Choose starting points for the orthogonal (to the contour) gridlines + # Use the existing cells of the standard mesh + all_cell_subset = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 5) + all_border_edges = SOLPS2IMAS.get_subset_boundary(space, all_cell_subset) + core_edges = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 15) + outer_target = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 13) + inner_target = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 14) + ci = [core_edges.element[i].object[1].index for i ∈ 1:length(core_edges.element)] + oi = + [outer_target.element[i].object[1].index for i ∈ 1:length(outer_target.element)] + ii = + [inner_target.element[i].object[1].index for i ∈ 1:length(inner_target.element)] + border_edges = [] + for i ∈ 1:length(all_border_edges) + bi = all_border_edges[i].object[1].index + if !(bi in oi) & !(bi in ii) & !(bi in ci) + border_edges = [border_edges; all_border_edges[i]] + end + end + + npol = length(border_edges) + mesh_r = zeros((npol, nlvl)) + mesh_z = zeros((npol, nlvl)) + corner_idx = 1 + for i ∈ 1:npol + edge_idx = border_edges[i].object[1].index + node_idx = space.objects_per_dimension[2].object[edge_idx].nodes[corner_idx] + geo = space.objects_per_dimension[1].object[node_idx].geometry + mesh_r[i, 1] = geo[1] + mesh_z[i, 1] = geo[2] + end + + # Step along the paths of steepest descent to populate the mesh. + dpsindr, dpsindz = OMAS.gradient(r_eq, z_eq, psin_eq) + psin_int = Interpolations.linear_interpolation((r_eq, z_eq), psin_eq) + dpdr = Interpolations.linear_interpolation((r_eq, z_eq), dpsindr) + dpdz = Interpolations.linear_interpolation((r_eq, z_eq), dpsindz) + rlim = (minimum(r_eq), maximum(r_eq)) + zlim = (minimum(z_eq), maximum(z_eq)) + pfr = psin_int.(mesh_r[:, 1], mesh_z[:, 1]) .< 1 + for i ∈ 1:npol + if (mesh_r[i, 1] > rlim[1]) & (mesh_r[i, 1] < rlim[2]) & + (mesh_z[i, 1] > zlim[1]) & (mesh_z[i, 1] < zlim[2]) + # direction = psin_int(mesh_r[i, 1], mesh_z[i, 1]) >= 1 ? 1 : -1 + direction = pfr[i] ? -1 : 1 + for j ∈ 2:nlvl + # This is low resolution linear shooting that could go wrong + mesh_r[i, j] = mesh_r[i, j-1] + mesh_z[i, j] = mesh_z[i, j-1] + upscale = 5 # Should improve gradient following in regions of low gradient (near X-point) + for k ∈ 1:upscale + if (mesh_r[i, j] > rlim[1]) & (mesh_r[i, j] < rlim[2]) & + (mesh_z[i, j] > zlim[1]) & (mesh_z[i, j] < zlim[2]) + dpr = dpdr(mesh_r[i, j], mesh_z[i, j]) + dpz = dpdz(mesh_r[i, j], mesh_z[i, j]) + d = dpsin * direction / upscale / (dpr .^ 2 + dpz .^ 2) + mesh_r[i, j] += d * dpr + mesh_z[i, j] += d * dpz + end + end + end + end + # println("i=",i,"; r=",mesh_r[i, 1],":", mesh_r[i, end],",z=",mesh_z[i, 1], ":", mesh_z[i, end]) + end + + # There's a special path; the one that probably should've gone through the + # secondary X-point (if there is one). Any numerical error will make this + # path miss the X-point and go off to somewhere crazy instead, so instead of + # that, let's draw a straight line from the edge of the SOLPS mesh to the + # X-point. + bssx = eqt.boundary_secondary_separatrix.x_point + nsx = length(bssx) + println("there are ", nsx, " secondary x points") + if nsx > 0 + # There is hopefully only one X point on the secondary separatrix, but + # just in case some joker made a secondary snowflake or something crazy, + # we should pick which one we like best. + if nsx == 1 + xidx = 1 # Easy peasy + else + score = zeros(nsx) + rsx = [bssx[i].r for i ∈ 1:nsx] + zsx = [bssx[i].z for i ∈ 1:nsx] + # Being close to the vertically flipped position of the primary X-pt + # seems nice; let's reward that + bsx = eqt.boundary_separatrix.x_point + npx = length(bsx) + rpx = [bsx[i].r for i ∈ 1:npx] + zpx = [bsx[i].z for i ∈ 1:npx] + # Well, if we have a primary snowflake, there could be multiple + # primary X-points, so now we have to pick our favorite primary. + # Being close to the boundary is a good sign + min_ds2 = ones(npx) * 1e6 + for j ∈ 1:nsx + dr = [mesh_r[i, 1] - rpx[j] for i ∈ 1:npol] + dz = [mesh_z[i, 1] - zpx[j] for i ∈ 1:npol] + ds2 = dr .^ 2 .+ dz .^ 2 + min_ds2[j] = minimum(ds2) + end + primary_idx = argmin(min_ds2) + rpx = rpx[primary_idx] + zpx = zpx[primary_idx] + # There can be only one + + # Distance between secondary X-points and vert flip of favorite primary + ds2xx = (rsx .- rpx) .^ 2 .+ (zsx .- (-zpx)) .^ 2 + score += ds2xx * 5 + + # Being close to the boundary is good + min_ds2xb = ones(nsx) * 1e6 + for j ∈ 1:nsx + dr = [mesh_r[i, 1] - rsx[j] for i ∈ 1:npol] + dz = [mesh_z[i, 1] - zsx[j] for i ∈ 1:npol] + ds2 = dr .^ 2 .+ dz .^ 2 + min_ds2xb[j] = minimum(ds2) + end + score += min_ds2xb + # And the winner is the one with the lowest score + xidx = argmin(score) + end + # Pick which edges are closest to the X-point + rx = bssx[xidx].r + zx = bssx[xidx].z + ds2 = (mesh_r[:, 1] .- rx) .^ 2 .+ (mesh_z[:, 1] .- zx) .^ 2 + closest = argmin(ds2) + # Work out the mesh spacing + drm = mesh_r[closest, 2] - mesh_r[closest, 1] + dzm = mesh_z[closest, 2] - mesh_z[closest, 1] + dm = sqrt(drm^2 + dzm^2) + # Pick the angle + dr = rx - mesh_r[closest, 1] + dz = zx - mesh_z[closest, 1] + angle = atan(dz, dr) + # Replace the points + for i ∈ 2:nlvl + mesh_r[closest, i] = mesh_r[closest, 1] + dm * (i - 1) * cos(angle) + mesh_z[closest, i] = mesh_z[closest, 1] + dm * (i - 1) * sin(angle) + end + end + + # ---------------------------------------------------------------------------------- + # Now we have all the nodes needed for the new mesh, but none are connected + # yet. They are, however, organized nicely in order, so it shouldn't be too + # hard. Any X-points in the domain should be aggressively ignored, so all the + # connections should be nice and simple. + # The PFR flag (that was calculated earlier) needs to be respected; pfr nodes + # don't connect to non-pfr nodes + o0 = space.objects_per_dimension[1] # Nodes + o1 = space.objects_per_dimension[2] # Edges + o2 = space.objects_per_dimension[3] # Cells (2D projection) + + n_old_node = length(o0.object) + n_new_node = nlvl * npol + n_nodes = n_old_node + n_new_node + + n_old_edge = length(o1.object) + n_new_edge_i = nlvl * (npol - 2) + n_new_edge_j = (nlvl - 1) * npol + n_new_edge = n_new_edge_i + n_new_edge_j + n_edges = n_old_edge + n_new_edge + + n_old_cell = length(o2.object) + n_new_cell = (nlvl - 1) * (npol - 2) # -2 to account for disconnect @ pfr + n_cells = n_old_cell + n_new_cell + + pfr_transition = argmax(abs.(diff(pfr))) + + # Define starting points and increments + node_start = n_old_node + 1 + edge_start1 = n_old_edge + 1 + edge_start2 = edge_start1 + n_new_edge_i + cell_start = n_old_cell + 1 + + n_per_i = nlvl + n_per_j = 1 + e1_per_i = nlvl # Edges along the npol direction / i direction + e2_per_i = nlvl - 1 # Edges along the nlvl direction / j direction + c_per_i = nlvl - 1 # # of cells < # of nodes + + # Get new subsets ready + n_existing_subsets = length(grid_ggd.grid_subset) + new_subset_indices = [-1, -2, -3, -4, -5, -201, -202, -203, -204, -205] + n_new_subsets = length(new_subset_indices) + resize!(grid_ggd.grid_subset, n_existing_subsets + n_new_subsets) + for i ∈ 1:n_new_subsets + grid_ggd.grid_subset[n_existing_subsets+i].identifier.index = + new_subset_indices[i] + end + ext_nodes_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -201) + ext_edges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -202) + ext_xedges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -203) + ext_yedges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -204) + ext_cells_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -205) + # resize!(ext_nodes_sub.element, n_new_node) + # resize!(ext_edges_sub.element, n_new_edge) + # resize!(ext_xedges_sub.element, n_new_edge_i) + # resize!(ext_yedges_sub.element, n_new_edge_j) + # resize!(ext_cells_sub.element, n_new_cell) + + # Preserve record of standard (non extended) mesh + for i ∈ 1:5 + std_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -i) + orig_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, i) + resize!(std_sub.element, length(orig_sub.element)) + for j ∈ 1:length(orig_sub.element) + std_sub.element[j] = deepcopy(orig_sub.element[j]) + end + std_sub.identifier.index = -i + std_sub.dimension = deepcopy(orig_sub.dimension) + std_sub.metric = deepcopy(orig_sub.metric) + end + all_nodes_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 1) + all_edges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 2) + all_xedges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 3) + all_yedges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 4) + all_cells_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 5) + # resize!(all_nodes_sub.element, n_nodes) + # resize!(all_edges_sub.element, n_edges) + # resize!(all_xedges_sub.element, length(all_xedges_sub.element) + n_new_edge_i) + # resize!(all_yedges_sub.element, length(all_yedges_sub.element) + n_new_edge_j) + # resize!(all_cells_sub.element, n_cells) + + nodes = resize!(o0.object, n_nodes) + edges = resize!(o1.object, n_edges) + cells = resize!(o2.object, n_cells) + for i ∈ 1:npol + pastpc = i > pfr_transition + for j ∈ 1:nlvl + # Modified counters + ii = i - 1 # Offset due to indexing from 1 + jj = j - 1 # Offset due to indexing from 1 + iii = ii - 1 - pastpc # # of connections < than # of nodes, missing row @ PFR transition + jjj = jj - 1 # # of connections < # of nodes + + # Nodes + node_idx = node_start + ii * n_per_i + jj + nodes[node_idx].geometry = [mesh_r[i, j], mesh_z[i, j]] + SOLPS2IMAS.add_subset_element!(ext_nodes_sub, space_idx, 0, node_idx) + SOLPS2IMAS.add_subset_element!(all_nodes_sub, space_idx, 0, node_idx) + + # Edges + if (i > 1) & (i != pfr_transition) # i-1 to i in the npol direction + edge_idx1 = edge_start1 + iii * e1_per_i + jj + edges[edge_idx1].nodes = [node_idx, node_idx - n_per_i] + SOLPS2IMAS.add_subset_element!(ext_edges_sub, space_idx, 1, edge_idx1) + SOLPS2IMAS.add_subset_element!(ext_xedges_sub, space_idx, 1, edge_idx1) + SOLPS2IMAS.add_subset_element!(all_edges_sub, space_idx, 1, edge_idx1) + SOLPS2IMAS.add_subset_element!(all_xedges_sub, space_idx, 1, edge_idx1) + end + if (j > 1) # j-1 to j in the nlvl direction + edge_idx2 = edge_start2 + ii * e2_per_i + jjj + edges[edge_idx2].nodes = [node_idx, node_idx - n_per_j] + SOLPS2IMAS.add_subset_element!(ext_edges_sub, space_idx, 1, edge_idx2) + SOLPS2IMAS.add_subset_element!(ext_yedges_sub, space_idx, 1, edge_idx2) + SOLPS2IMAS.add_subset_element!(all_edges_sub, space_idx, 1, edge_idx2) + SOLPS2IMAS.add_subset_element!(all_yedges_sub, space_idx, 1, edge_idx2) + end + + # Cells + if (i > 1) & (i != pfr_transition) & (j > 1) + cell_idx = cell_start + iii * c_per_i + jjj + cells[cell_idx].nodes = [ + node_idx, + node_idx - n_per_i, + node_idx - n_per_i - n_per_j, + node_idx - n_per_j, + ] + SOLPS2IMAS.add_subset_element!(ext_cells_sub, space_idx, 2, cell_idx) + SOLPS2IMAS.add_subset_element!(all_cells_sub, space_idx, 2, cell_idx) + end + end + end + + # # Make the contours + # for psin_level ∈ psin_levels + # if psin_level >= secondary_psin + # throw( + # ArgumentError( + # "a contour level is outside the secondary separatrix, which isn't handled yet", + # ), + # ) + # end + # c = Contour.contour(r_eq, z_eq, psin_eq, psin_level) + # clines = Contour.lines(c) + # end + fr = open("mesh_r.dat", "w") + fz = open("mesh_z.dat", "w") + for i ∈ 1:npol + for j ∈ 1:nlvl + print(fr, mesh_r[i, j], " ") + print(fz, mesh_z[i, j], " ") + end + println(fr, "") + println(fz, "") + end + return +end + """ fill_in_extrapolated_edge_profile!( dd::OMAS.dd, quantity_name::String; method::String="simple", diff --git a/test/runtests.jl b/test/runtests.jl index 116fcf8..52bb24a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -108,7 +108,7 @@ function define_default_sample_set() b2fgmtry, b2time, b2mn, gridspec, eqdsk = file_list eqdsk = splitdir(pathof(SD4SOLPS))[1] * - "/../sample/ITER_Lore_2296_00000/EQDSK/Baseline2008-li0.70.x4.mod2.eqdsk" + "/../sample/ITER_Lore_2296_00000/EQDSK/g002296.00200" return b2fgmtry, b2time, b2mn, gridspec, eqdsk end @@ -225,6 +225,43 @@ if args["edge_profile_extension"] dpsin = SD4SOLPS.mesh_psi_spacing(dd) @test dpsin > 0.0 + # Extend the mesh + grid_ggd_idx = 1 + grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] + extended_subs = 1:5 + orig_subs = [ + deepcopy(SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, i)) for + i ∈ extended_subs + ] + SD4SOLPS.mesh_extension_sol!(dd; grid_ggd_idx=grid_ggd_idx) + for j ∈ extended_subs + orig_sub = orig_subs[j] + std_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -j) + all_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, j) + ext_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -200 - j) + no = length(orig_sub.element) + ns = length(std_sub.element) + na = length(all_sub.element) + ne = length(ext_sub.element) + orig_indices = [orig_sub.element[i].object[1].index for i ∈ 1:no] + std_indices = [std_sub.element[i].object[1].index for i ∈ 1:ns] + all_indices = [all_sub.element[i].object[1].index for i ∈ 1:na] + ext_indices = [ext_sub.element[i].object[1].index for i ∈ 1:ne] + @test std_sub.identifier.index == -j + @test all(orig_indices .== std_indices) + all_indices_reconstruct = [std_indices; ext_indices] + @test all(all_indices .== all_indices_reconstruct) + + # Verify that original and standard are separate and not refs to each other + arbitrary_change = 5 + arb_el = 2 + orig_sub.element[arb_el].object[1].index += arbitrary_change + @test orig_sub.element[arb_el].object[1].index != + std_sub.element[arb_el].object[1].index + orig_sub.element[arb_el].object[1].index -= arbitrary_change + end + + # Prepare profiles that need to be extended n_edge = 47 n_outer_prof = 13 quantity_edge = From 131e20402ea3b4eabb3af3a1bd61406a0ccfa7b9 Mon Sep 17 00:00:00 2001 From: David Eldon Date: Fri, 3 Nov 2023 13:43:28 -0700 Subject: [PATCH 05/12] Break mesh extension function into smaller functions --- src/supersize_profile.jl | 191 +++++++++++++++++++++++---------------- 1 file changed, 115 insertions(+), 76 deletions(-) diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index 333f4bc..3c46f2f 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -280,6 +280,19 @@ function extrapolate_edge_exp( return y0 * exp(-x ./ lambda) end +function prep_flux_map(dd::OMAS.dd; eq_time_idx::Int64=1, eq_profiles_2d_idx::Int64=1) + eqt = dd.equilibrium.time_slice[eq_time_idx] + p2 = eqt.profiles_2d[eq_profiles_2d_idx] + r_eq = p2.grid.dim1 + z_eq = p2.grid.dim2 + psi = p2.psi + psia = eqt.global_quantities.psi_axis + psib = eqt.global_quantities.psi_boundary + psin_eq = (psi .- psia) ./ (psib - psia) + rzpi = Interpolations.LinearInterpolation((r_eq, z_eq), psin_eq) + return r_eq, z_eq, psin_eq, rzpi +end + #! format off """ mesh_psi_spacing( @@ -342,17 +355,7 @@ function mesh_psi_spacing( end # Get flux map - eqt = dd.equilibrium.time_slice[eq_time_idx] - p2 = eqt.profiles_2d[eq_profiles_2d_idx] - r_eq = p2.grid.dim1 - z_eq = p2.grid.dim2 - psi = p2.psi - psia = eqt.global_quantities.psi_axis - psib = eqt.global_quantities.psi_boundary - psin_eq = (psi .- psia) ./ (psib - psia) - rzpi = Interpolations.linear_interpolation((r_eq, z_eq), psin_eq) - # println(minimum(r_eq), ", ", maximum(r_eq)) - # println(minimum(z_eq), ", ", maximum(z_eq)) + r_eq, z_eq, psin_eq, rzpi = prep_flux_map(dd; eq_time_idx, eq_profiles_2d_idx) # Get a row of cells. Since the mesh should be aligned to the flux surfaces, # it shouldn't matter which row is used, although the divertor rows might be @@ -385,28 +388,14 @@ function mesh_psi_spacing( return dpsin end -""" - function mesh_extension_sol() - -Extends the mesh out into the SOL -""" -function mesh_extension_sol!( +function pick_extension_psi_range( dd::OMAS.dd; eq_time_idx::Int64=1, eq_profiles_2d_idx::Int64=1, grid_ggd_idx::Int64=1, space_idx::Int64=1, ) - # Prep flux map - eqt = dd.equilibrium.time_slice[eq_time_idx] - p2 = eqt.profiles_2d[eq_profiles_2d_idx] - r_eq = p2.grid.dim1 - z_eq = p2.grid.dim2 - psi = p2.psi - psia = eqt.global_quantities.psi_axis - psib = eqt.global_quantities.psi_boundary - psin_eq = (psi .- psia) ./ (psib - psia) - rzpi = Interpolations.LinearInterpolation((r_eq, z_eq), psin_eq) + r_eq, z_eq, psin_eq, rzpi = prep_flux_map(dd; eq_time_idx, eq_profiles_2d_idx) # Use wall to find maximum extent of contouring project limiter = dd.wall.description_2d[1].limiter @@ -417,8 +406,7 @@ function mesh_extension_sol!( # Use ggd mesh to find inner limit of contouring project grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] space = grid_ggd.space[space_idx] - midplane_subset = - SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 11) + midplane_subset = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 11) midplane_cell_centers = GGDUtils.get_subset_centers(space, midplane_subset) psin_midplane = rzpi.(midplane_cell_centers[end][1], midplane_cell_centers[end][2]) @@ -438,12 +426,23 @@ function mesh_extension_sol!( lvlend = maximum(wall_psin * sign(dpsin)) / sign(dpsin) + dpsin nlvl = Int64(ceil((lvlend - lvlstart) / dpsin)) psin_levels = collect(LinRange(lvlstart, lvlend, nlvl)) - if hasproperty(eqt.boundary_secondary_separatrix, :psi) - secondary_psi = eqt.boundary_secondary_separatrix.psi - secondary_psin = (secondary_psi - psia) / (psib - psia) - else - secondary_psin = lvlend + dpsin - end + # eqt = dd.equilibrium.time_slice[eq_time_idx] + # if hasproperty(eqt.boundary_secondary_separatrix, :psi) + # secondary_psi = eqt.boundary_secondary_separatrix.psi + # secondary_psin = (secondary_psi - psia) / (psib - psia) + # else + # secondary_psin = lvlend + dpsin + # end + return psin_levels +end + +function pick_mesh_ext_starting_points( + dd::OMAS.dd; + grid_ggd_idx::Int64=1, + space_idx::Int64=1, +) + grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] + space = grid_ggd.space[space_idx] # Choose starting points for the orthogonal (to the contour) gridlines # Use the existing cells of the standard mesh @@ -466,29 +465,39 @@ function mesh_extension_sol!( end npol = length(border_edges) - mesh_r = zeros((npol, nlvl)) - mesh_z = zeros((npol, nlvl)) + r = zeros(npol) + z = zeros(npol) corner_idx = 1 for i ∈ 1:npol edge_idx = border_edges[i].object[1].index node_idx = space.objects_per_dimension[2].object[edge_idx].nodes[corner_idx] geo = space.objects_per_dimension[1].object[node_idx].geometry - mesh_r[i, 1] = geo[1] - mesh_z[i, 1] = geo[2] + r[i, 1] = geo[1] + z[i, 1] = geo[2] + end + return r, z +end + +function mesh_ext_follow_grad(r_eq, z_eq, psin_eq, rzpi, rstart, zstart, nlvl, dpsin) + npol = length(rstart) + mesh_r = zeros((npol, nlvl)) + mesh_z = zeros((npol, nlvl)) + for i ∈ 1:npol + mesh_r[i, 1] = rstart[i] + mesh_z[i, 1] = zstart[i] end # Step along the paths of steepest descent to populate the mesh. dpsindr, dpsindz = OMAS.gradient(r_eq, z_eq, psin_eq) - psin_int = Interpolations.linear_interpolation((r_eq, z_eq), psin_eq) dpdr = Interpolations.linear_interpolation((r_eq, z_eq), dpsindr) dpdz = Interpolations.linear_interpolation((r_eq, z_eq), dpsindz) rlim = (minimum(r_eq), maximum(r_eq)) zlim = (minimum(z_eq), maximum(z_eq)) - pfr = psin_int.(mesh_r[:, 1], mesh_z[:, 1]) .< 1 + pfr = rzpi.(mesh_r[:, 1], mesh_z[:, 1]) .< 1 for i ∈ 1:npol if (mesh_r[i, 1] > rlim[1]) & (mesh_r[i, 1] < rlim[2]) & (mesh_z[i, 1] > zlim[1]) & (mesh_z[i, 1] < zlim[2]) - # direction = psin_int(mesh_r[i, 1], mesh_z[i, 1]) >= 1 ? 1 : -1 + # direction = rzpi(mesh_r[i, 1], mesh_z[i, 1]) >= 1 ? 1 : -1 direction = pfr[i] ? -1 : 1 for j ∈ 2:nlvl # This is low resolution linear shooting that could go wrong @@ -509,12 +518,16 @@ function mesh_extension_sol!( end # println("i=",i,"; r=",mesh_r[i, 1],":", mesh_r[i, end],",z=",mesh_z[i, 1], ":", mesh_z[i, end]) end + return mesh_r, mesh_z +end +function modify_mesh_ext_near_x!(eqt, mesh_r, mesh_z) # There's a special path; the one that probably should've gone through the # secondary X-point (if there is one). Any numerical error will make this # path miss the X-point and go off to somewhere crazy instead, so instead of # that, let's draw a straight line from the edge of the SOLPS mesh to the # X-point. + npol, nlvl = size(mesh_r) bssx = eqt.boundary_secondary_separatrix.x_point nsx = length(bssx) println("there are ", nsx, " secondary x points") @@ -584,14 +597,13 @@ function mesh_extension_sol!( mesh_z[closest, i] = mesh_z[closest, 1] + dm * (i - 1) * sin(angle) end end +end - # ---------------------------------------------------------------------------------- - # Now we have all the nodes needed for the new mesh, but none are connected - # yet. They are, however, organized nicely in order, so it shouldn't be too - # hard. Any X-points in the domain should be aggressively ignored, so all the - # connections should be nice and simple. - # The PFR flag (that was calculated earlier) needs to be respected; pfr nodes - # don't connect to non-pfr nodes +function record_regular_mesh!(dd::OMAS.dd, grid_ggd_idx, space_idx, mesh_r, mesh_z, cut) + grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] + space = grid_ggd.space[space_idx] + + npol, nlvl = size(mesh_r) o0 = space.objects_per_dimension[1] # Nodes o1 = space.objects_per_dimension[2] # Edges o2 = space.objects_per_dimension[3] # Cells (2D projection) @@ -610,8 +622,6 @@ function mesh_extension_sol!( n_new_cell = (nlvl - 1) * (npol - 2) # -2 to account for disconnect @ pfr n_cells = n_old_cell + n_new_cell - pfr_transition = argmax(abs.(diff(pfr))) - # Define starting points and increments node_start = n_old_node + 1 edge_start1 = n_old_edge + 1 @@ -638,11 +648,6 @@ function mesh_extension_sol!( ext_xedges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -203) ext_yedges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -204) ext_cells_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -205) - # resize!(ext_nodes_sub.element, n_new_node) - # resize!(ext_edges_sub.element, n_new_edge) - # resize!(ext_xedges_sub.element, n_new_edge_i) - # resize!(ext_yedges_sub.element, n_new_edge_j) - # resize!(ext_cells_sub.element, n_new_cell) # Preserve record of standard (non extended) mesh for i ∈ 1:5 @@ -661,17 +666,12 @@ function mesh_extension_sol!( all_xedges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 3) all_yedges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 4) all_cells_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 5) - # resize!(all_nodes_sub.element, n_nodes) - # resize!(all_edges_sub.element, n_edges) - # resize!(all_xedges_sub.element, length(all_xedges_sub.element) + n_new_edge_i) - # resize!(all_yedges_sub.element, length(all_yedges_sub.element) + n_new_edge_j) - # resize!(all_cells_sub.element, n_cells) nodes = resize!(o0.object, n_nodes) edges = resize!(o1.object, n_edges) cells = resize!(o2.object, n_cells) for i ∈ 1:npol - pastpc = i > pfr_transition + pastpc = i > cut for j ∈ 1:nlvl # Modified counters ii = i - 1 # Offset due to indexing from 1 @@ -686,7 +686,7 @@ function mesh_extension_sol!( SOLPS2IMAS.add_subset_element!(all_nodes_sub, space_idx, 0, node_idx) # Edges - if (i > 1) & (i != pfr_transition) # i-1 to i in the npol direction + if (i > 1) & (i != cut) # i-1 to i in the npol direction edge_idx1 = edge_start1 + iii * e1_per_i + jj edges[edge_idx1].nodes = [node_idx, node_idx - n_per_i] SOLPS2IMAS.add_subset_element!(ext_edges_sub, space_idx, 1, edge_idx1) @@ -704,7 +704,7 @@ function mesh_extension_sol!( end # Cells - if (i > 1) & (i != pfr_transition) & (j > 1) + if (i > 1) & (i != cut) & (j > 1) cell_idx = cell_start + iii * c_per_i + jjj cells[cell_idx].nodes = [ node_idx, @@ -718,18 +718,6 @@ function mesh_extension_sol!( end end - # # Make the contours - # for psin_level ∈ psin_levels - # if psin_level >= secondary_psin - # throw( - # ArgumentError( - # "a contour level is outside the secondary separatrix, which isn't handled yet", - # ), - # ) - # end - # c = Contour.contour(r_eq, z_eq, psin_eq, psin_level) - # clines = Contour.lines(c) - # end fr = open("mesh_r.dat", "w") fz = open("mesh_z.dat", "w") for i ∈ 1:npol @@ -740,6 +728,57 @@ function mesh_extension_sol!( println(fr, "") println(fz, "") end +end + +""" + function mesh_extension_sol() + +Extends the mesh out into the SOL +""" +function mesh_extension_sol!( + dd::OMAS.dd; + eq_time_idx::Int64=1, + eq_profiles_2d_idx::Int64=1, + grid_ggd_idx::Int64=1, + space_idx::Int64=1, +) + grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] + space = grid_ggd.space[space_idx] + eqt = dd.equilibrium.time_slice[eq_time_idx] + + r_eq, z_eq, psin_eq, rzpi = prep_flux_map(dd; eq_time_idx, eq_profiles_2d_idx) + psin_levels = pick_extension_psi_range( + dd; + eq_time_idx, + eq_profiles_2d_idx, + grid_ggd_idx, + space_idx, + ) + nlvl = length(psin_levels) + dpsin = psin_levels[2] - psin_levels[1] + grad_start_r, grad_start_z = + pick_mesh_ext_starting_points(dd; grid_ggd_idx, space_idx) + npol = length(grad_start_r) + mesh_r, mesh_z = mesh_ext_follow_grad( + r_eq, + z_eq, + psin_eq, + rzpi, + grad_start_r, + grad_start_z, + nlvl, + dpsin, + ) + modify_mesh_ext_near_x!(eqt, mesh_r, mesh_z) + + # Now we have all the nodes needed for the new mesh, but none are connected + # yet. They are, however, organized nicely in order, so it shouldn't be too + # hard. Any X-points in the domain should be aggressively ignored, so all the + # connections should be nice and simple. + # The PFR flag needs to be respected; pfr nodes don't connect to non-pfr nodes + pfr = rzpi.(mesh_r[:, 1], mesh_z[:, 1]) .< 1 + pfr_transition = argmax(abs.(diff(pfr))) + record_regular_mesh!(dd, grid_ggd_idx, space_idx, mesh_r, mesh_z, pfr_transition) return end From 8b60028646ade392b8fe92cd50070fd8de4977d3 Mon Sep 17 00:00:00 2001 From: David Eldon Date: Mon, 6 Nov 2023 09:32:03 -0800 Subject: [PATCH 06/12] Add caching for mesh extensions - Add documentation to edge mesh extension functions - JSON instead of YAML: turning the matrices back into matrices is easier this way --- Project.toml | 1 + data/.gitignore | 2 + src/SD4SOLPS.jl | 1 + src/supersize_profile.jl | 218 +++++++++++++++++++++++++++++++++++++-- test/runtests.jl | 4 +- 5 files changed, 214 insertions(+), 12 deletions(-) create mode 100644 data/.gitignore diff --git a/Project.toml b/Project.toml index e1d5140..f727ad2 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ Contour = "d38c429a-6771-53c6-b99e-75d170b6e991" EFIT = "cda752c5-6b03-55a3-9e33-132a441b0c17" GGDUtils = "b7b5e640-9b39-4803-84eb-376048795def" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" OMAS = "91cfaa06-6526-4804-8666-b540b3feef2f" PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab" PlotUtils = "995b91a9-d308-5afd-9ec6-746e21dbc043" diff --git a/data/.gitignore b/data/.gitignore new file mode 100644 index 0000000..f6bc585 --- /dev/null +++ b/data/.gitignore @@ -0,0 +1,2 @@ +*.mesh_ext.json + diff --git a/src/SD4SOLPS.jl b/src/SD4SOLPS.jl index 09a9f8e..0f87fe4 100644 --- a/src/SD4SOLPS.jl +++ b/src/SD4SOLPS.jl @@ -280,6 +280,7 @@ function preparation( # ... more profiles here as they become available in b2time println("Extrapolated core profiles") + cached_mesh_extension!(dd, eqdsk_file, b2fgmtry) fill_in_extrapolated_edge_profile!(dd, "electrons.density"; method=core_method) # ... more profiles here println("Extrapolated edge profiles (but not really (placeholder only))") diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index 3c46f2f..33d19af 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -8,6 +8,7 @@ using Interpolations: Interpolations using GGDUtils: GGDUtils using PolygonOps: PolygonOps using SOLPS2IMAS: SOLPS2IMAS +using JSON: JSON export extrapolate_core export fill_in_extrapolated_core_profile @@ -280,6 +281,18 @@ function extrapolate_edge_exp( return y0 * exp(-x ./ lambda) end +""" + prep_flux_map() + +Reads equilibrium data and extracts/derives some useful quantities. +This is very basic, but it was being repeated and that's a no-no. +Returns: + + - R values of the equilibrium grid + - Z values of the eq grid + - normalized poloidal flux on the equilibrium grid + - a linear interpolation of norm pol flux vs. R and Z, ready to be evaluated +""" function prep_flux_map(dd::OMAS.dd; eq_time_idx::Int64=1, eq_profiles_2d_idx::Int64=1) eqt = dd.equilibrium.time_slice[eq_time_idx] p2 = eqt.profiles_2d[eq_profiles_2d_idx] @@ -388,6 +401,14 @@ function mesh_psi_spacing( return dpsin end +""" + pick_extension_psi_range() + +Defines the psi_N levels for an extended mesh. The range of psi_N levels starts +at the outer edge of the existing edge_profiles mesh at the midplane and goes +out to the most distant (in flux space) point on the limiting surface. +Returns a vector of psi_N levels. +""" function pick_extension_psi_range( dd::OMAS.dd; eq_time_idx::Int64=1, @@ -436,6 +457,17 @@ function pick_extension_psi_range( return psin_levels end +""" + pick_mesh_ext_starting_points(dd; grid_ggd_idx, space_idx) + +Picks starting points for the radial lines of the mesh extension. The strategy +is to start from the outer edge of the existing mesh and follow the steepest +gradient (of psi_N) to extend these gridlines outward. +dd: a data dictionary instance with edge_profiles ggd and equilibrium loaded +grid_ggd_idx: index within ggd +space_idx: space number / index of the space to work with within edge_profiles +Returns a tuple with vectors of R and Z starting points. +""" function pick_mesh_ext_starting_points( dd::OMAS.dd; grid_ggd_idx::Int64=1, @@ -478,7 +510,38 @@ function pick_mesh_ext_starting_points( return r, z end -function mesh_ext_follow_grad(r_eq, z_eq, psin_eq, rzpi, rstart, zstart, nlvl, dpsin) +#!format off +""" + mesh_ext_follow_grad() + +Follows the steepest gradient from a set of starting points, dropping nodes at +approximately regular intervals in psi_N. Due to the numerical techniques used, the +node spacing may be imperfect (especially prone to error in regions where curvature +of psi_N is large compared to its gradient). +r_eq: Equilibrium reconstruction's grid, R coordinates +z_eq: Equilibrium reconstruction's grid, Z coordinates +psin_eq: Normalized poloidal flux in the equilibrium reconstruction as a function of R and Z +rzpi: linear interpolation of psin_eq() as a function of r_eq and z_eq + This was probably already computed and I think time would be saved by reusing it. + If you don't already have it, you can pass in nothing and let this function calculate it. +rstart: R coordinates of starting points for the gradient following. +zstart: Z coordinates of starting points +nlvl: number of nodes to drop while following the gradient +dpsin: node spacing in delta psi_N + +Returns two matrices with R and Z coordinates of the mesh extension +""" +#!format on +function mesh_ext_follow_grad( + r_eq::Vector{Float64}, + z_eq::Vector{Float64}, + psin_eq::Matrix, + rzpi, + rstart::Vector{Float64}, + zstart::Vector{Float64}, + nlvl::Int64, + dpsin::Float64, +) npol = length(rstart) mesh_r = zeros((npol, nlvl)) mesh_z = zeros((npol, nlvl)) @@ -487,6 +550,10 @@ function mesh_ext_follow_grad(r_eq, z_eq, psin_eq, rzpi, rstart, zstart, nlvl, d mesh_z[i, 1] = zstart[i] end + if rzpi === nothing + rzpi = Interpolations.LinearInterpolation((r_eq, z_eq), psin_eq) + end + # Step along the paths of steepest descent to populate the mesh. dpsindr, dpsindz = OMAS.gradient(r_eq, z_eq, psin_eq) dpdr = Interpolations.linear_interpolation((r_eq, z_eq), dpsindr) @@ -521,6 +588,15 @@ function mesh_ext_follow_grad(r_eq, z_eq, psin_eq, rzpi, rstart, zstart, nlvl, d return mesh_r, mesh_z end +""" + modify_mesh_ext_near_x! + +Modifies an extended mesh near a secondary X-point to compensate for the +tendency of the mesh to go nuts near the X-point. +eqt: equilibrium.time_slice information +mesh_r: matrix of R values for the extended mesh +mesh_z: matrix of Z values for the extended mesh +""" function modify_mesh_ext_near_x!(eqt, mesh_r, mesh_z) # There's a special path; the one that probably should've gone through the # secondary X-point (if there is one). Any numerical error will make this @@ -599,7 +675,30 @@ function modify_mesh_ext_near_x!(eqt, mesh_r, mesh_z) end end -function record_regular_mesh!(dd::OMAS.dd, grid_ggd_idx, space_idx, mesh_r, mesh_z, cut) +#!format off +""" + record_regular_mesh!() + +Records arrays of mesh data from regular 2D arrays into the DD +dd: the data dictionary +grid_ggd_idx: index of the grid_ggd within edge_profiles +space_idx: index of the space in edge_profiles +mesh_r: Matrix of R values along a mesh. Should be 2D. The two dimensions are in + the radial and poloidal directions. +mesh_z: Z values to go with mesh_r. +cut: Poloidal index of a cut between two groups of disconnected cells. Poloidal + connections (faces, cells) will not be added between this poloidal index + and the next index. +""" +#!format on +function record_regular_mesh!( + dd::OMAS.dd, + grid_ggd_idx::Int64, + space_idx::Int64, + mesh_r::Matrix{Float64}, + mesh_z::Matrix{Float64}, + cut::Int64, +) grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] space = grid_ggd.space[space_idx] @@ -717,17 +816,114 @@ function record_regular_mesh!(dd::OMAS.dd, grid_ggd_idx, space_idx, mesh_r, mesh end end end +end - fr = open("mesh_r.dat", "w") - fz = open("mesh_z.dat", "w") - for i ∈ 1:npol - for j ∈ 1:nlvl - print(fr, mesh_r[i, j], " ") - print(fz, mesh_z[i, j], " ") +""" + convert_filename(filename::String) + +Converts a filename into a string that doesn't have illegal characters. +The main application is removing the path separator from source files with full +paths so the full paths* can be part of the new filename. This way, the input +files used to form some data can be part of the cache name, allowing quick lookup: +the cache filename is defined by the input files, and if it doesn't exist, it +needs to be generated. +""" +function convert_filename(filename::String) + filename_mod = replace(filename, "/" => "__") # Illegal on *nix bc it's the path separator + filename_mod = replace(filename_mod, ":" => "--") # Illegal on mac and windows + filename_mod = replace(filename_mod, "\\" => "__") # Illegal on windows + filename_mod = replace(filename_mod, "\"" => "''") # Illegal on windows + filename_mod = replace(filename_mod, "<" => "_lt_") # Illegal on windows + filename_mod = replace(filename_mod, ">" => "_gt_") # Illegal on windows + filename_mod = replace(filename_mod, "|" => "_pipe_") # Illegal on windows + filename_mod = replace(filename_mod, "?" => "_q_") # Illegal on windows + filename_mod = replace(filename_mod, "*" => "_a_") # Illegal on windows + return filename_mod +end + +#!format off +""" + cached_mesh_extension!() + +Adds an extended mesh to a data dictionary, possibly from a cached result. +dd: The data dictionary. It will be modified in place. +eqdsk_file: the name of the EQDSK file that was used to get equilibrium data in + the dd. +b2fgmtry: the name of the SOLPS geometry file that was used to get GGD info in + edge_profiles in the dd. +eq_time_idx: Index of the time slice in equilibrium +eq_profiles_2d_idx: Index of the 2D profile set in equilibrium + (there is usually only one) +grid_ggd_idx: Index of the grid_ggd set in edge_profiles +space_idx: Index of the space +clear_cache: delete any existing cache file (for use in testing) +""" +#!format on +function cached_mesh_extension!( + dd::OMAS.dd, + eqdsk_file::String, + b2fgmtry::String; + eq_time_idx::Int64=1, + eq_profiles_2d_idx::Int64=1, + grid_ggd_idx::Int64=1, + space_idx::Int64=1, + clear_cache=false, +) + path = "$(@__DIR__)/../data/" + eqdsk_file_mod = convert_filename(eqdsk_file) + b2fgmtry_mod = convert_filename(eqdsk_file) + cached_ext_name = path * eqdsk_file_mod * "_" * b2fgmtry_mod * ".mesh_ext.json" + if clear_cache + rm(cached_ext_name; force=true) + return cached_ext_name + end + if isfile(cached_ext_name) + # md = YAML.load_file(cached_ext_name) + md = JSON.parsefile(cached_ext_name) + pfr_transition = md["pfr_transition"] + mesh_r = convert(Matrix{Float64}, mapreduce(permutedims, vcat, md["r"])') + mesh_z = convert(Matrix{Float64}, mapreduce(permutedims, vcat, md["z"])') + npol = md["npol"] + nlvl = md["nlvl"] + record_regular_mesh!( + dd, + grid_ggd_idx, + space_idx, + mesh_r, + mesh_z, + pfr_transition, + ) + else + mesh_r, mesh_z, pfr_transition = mesh_extension_sol!( + dd; + eq_time_idx=eq_time_idx, + eq_profiles_2d_idx=eq_profiles_2d_idx, + grid_ggd_idx=grid_ggd_idx, + space_idx=space_idx, + ) + npol, nlvl = size(mesh_r) + data = Dict() + data["pfr_transition"] = pfr_transition + data["npol"] = npol + data["nlvl"] = nlvl + data["r"] = mesh_r + data["z"] = mesh_z + # YAML.write_file(cached_ext_name, data) + open(cached_ext_name, "w") do f + JSON.print(f, data) end - println(fr, "") - println(fz, "") + # fr = open("mesh_r.dat", "w") + # fz = open("mesh_z.dat", "w") + # for i ∈ 1:npol + # for j ∈ 1:nlvl + # print(fr, mesh_r[i, j], " ") + # print(fz, mesh_z[i, j], " ") + # end + # println(fr, "") + # println(fz, "") + # end end + return cached_ext_name end """ @@ -779,7 +975,7 @@ function mesh_extension_sol!( pfr = rzpi.(mesh_r[:, 1], mesh_z[:, 1]) .< 1 pfr_transition = argmax(abs.(diff(pfr))) record_regular_mesh!(dd, grid_ggd_idx, space_idx, mesh_r, mesh_z, pfr_transition) - return + return mesh_r, mesh_z, pfr_transition end """ diff --git a/test/runtests.jl b/test/runtests.jl index 52bb24a..0547680 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -233,7 +233,9 @@ if args["edge_profile_extension"] deepcopy(SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, i)) for i ∈ extended_subs ] - SD4SOLPS.mesh_extension_sol!(dd; grid_ggd_idx=grid_ggd_idx) + cfn = SD4SOLPS.cached_mesh_extension!(dd, eqdsk, b2fgmtry; clear_cache=true) + println("cleared ext mesh cache: ", cfn) + SD4SOLPS.cached_mesh_extension!(dd, eqdsk, b2fgmtry; grid_ggd_idx=grid_ggd_idx) for j ∈ extended_subs orig_sub = orig_subs[j] std_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -j) From b6e205f532c248ba14774b1c3a0cb5bc078cea18 Mon Sep 17 00:00:00 2001 From: David Eldon Date: Mon, 6 Nov 2023 12:42:27 -0800 Subject: [PATCH 07/12] Adapt to GGD-related functions moving from SOLPS2IMAS to GGDUtils - See https://github.com/ProjectTorreyPines/GGDUtils.jl/pull/25 --- src/supersize_profile.jl | 71 ++++++++++++++++++++-------------------- test/runtests.jl | 9 ++--- 2 files changed, 41 insertions(+), 39 deletions(-) diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index 33d19af..bb237a3 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -5,7 +5,8 @@ Utilities for extrapolating profiles # import CalculusWithJulia using OMAS: OMAS using Interpolations: Interpolations -using GGDUtils: GGDUtils +using GGDUtils: + GGDUtils, get_grid_subset_with_index, add_subset_element!, get_subset_boundary using PolygonOps: PolygonOps using SOLPS2IMAS: SOLPS2IMAS using JSON: JSON @@ -136,9 +137,9 @@ function fill_in_extrapolated_core_profile!( grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] space = grid_ggd.space[space_idx] cell_subset = - SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 5) + get_grid_subset_with_index(grid_ggd, 5) midplane_subset = - SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 11) + get_grid_subset_with_index(grid_ggd, 11) if length(midplane_subset.element) < 1 throw( @@ -376,7 +377,7 @@ function mesh_psi_spacing( grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] space = grid_ggd.space[space_idx] midplane_subset = - SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 11) + get_grid_subset_with_index(grid_ggd, 11) midplane_cell_centers = GGDUtils.get_subset_centers(space, midplane_subset) r_mesh = [midplane_cell_centers[i][1] for i ∈ eachindex(midplane_cell_centers)] z_mesh = [midplane_cell_centers[i][2] for i ∈ eachindex(midplane_cell_centers)] @@ -427,7 +428,7 @@ function pick_extension_psi_range( # Use ggd mesh to find inner limit of contouring project grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] space = grid_ggd.space[space_idx] - midplane_subset = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 11) + midplane_subset = get_grid_subset_with_index(grid_ggd, 11) midplane_cell_centers = GGDUtils.get_subset_centers(space, midplane_subset) psin_midplane = rzpi.(midplane_cell_centers[end][1], midplane_cell_centers[end][2]) @@ -478,11 +479,11 @@ function pick_mesh_ext_starting_points( # Choose starting points for the orthogonal (to the contour) gridlines # Use the existing cells of the standard mesh - all_cell_subset = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 5) - all_border_edges = SOLPS2IMAS.get_subset_boundary(space, all_cell_subset) - core_edges = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 15) - outer_target = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 13) - inner_target = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 14) + all_cell_subset = get_grid_subset_with_index(grid_ggd, 5) + all_border_edges = get_subset_boundary(space, all_cell_subset) + core_edges = get_grid_subset_with_index(grid_ggd, 15) + outer_target = get_grid_subset_with_index(grid_ggd, 13) + inner_target = get_grid_subset_with_index(grid_ggd, 14) ci = [core_edges.element[i].object[1].index for i ∈ 1:length(core_edges.element)] oi = [outer_target.element[i].object[1].index for i ∈ 1:length(outer_target.element)] @@ -742,16 +743,16 @@ function record_regular_mesh!( grid_ggd.grid_subset[n_existing_subsets+i].identifier.index = new_subset_indices[i] end - ext_nodes_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -201) - ext_edges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -202) - ext_xedges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -203) - ext_yedges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -204) - ext_cells_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -205) + ext_nodes_sub = get_grid_subset_with_index(grid_ggd, -201) + ext_edges_sub = get_grid_subset_with_index(grid_ggd, -202) + ext_xedges_sub = get_grid_subset_with_index(grid_ggd, -203) + ext_yedges_sub = get_grid_subset_with_index(grid_ggd, -204) + ext_cells_sub = get_grid_subset_with_index(grid_ggd, -205) # Preserve record of standard (non extended) mesh for i ∈ 1:5 - std_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -i) - orig_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, i) + std_sub = get_grid_subset_with_index(grid_ggd, -i) + orig_sub = get_grid_subset_with_index(grid_ggd, i) resize!(std_sub.element, length(orig_sub.element)) for j ∈ 1:length(orig_sub.element) std_sub.element[j] = deepcopy(orig_sub.element[j]) @@ -760,11 +761,11 @@ function record_regular_mesh!( std_sub.dimension = deepcopy(orig_sub.dimension) std_sub.metric = deepcopy(orig_sub.metric) end - all_nodes_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 1) - all_edges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 2) - all_xedges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 3) - all_yedges_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 4) - all_cells_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 5) + all_nodes_sub = get_grid_subset_with_index(grid_ggd, 1) + all_edges_sub = get_grid_subset_with_index(grid_ggd, 2) + all_xedges_sub = get_grid_subset_with_index(grid_ggd, 3) + all_yedges_sub = get_grid_subset_with_index(grid_ggd, 4) + all_cells_sub = get_grid_subset_with_index(grid_ggd, 5) nodes = resize!(o0.object, n_nodes) edges = resize!(o1.object, n_edges) @@ -781,25 +782,25 @@ function record_regular_mesh!( # Nodes node_idx = node_start + ii * n_per_i + jj nodes[node_idx].geometry = [mesh_r[i, j], mesh_z[i, j]] - SOLPS2IMAS.add_subset_element!(ext_nodes_sub, space_idx, 0, node_idx) - SOLPS2IMAS.add_subset_element!(all_nodes_sub, space_idx, 0, node_idx) + add_subset_element!(ext_nodes_sub, space_idx, 0, node_idx) + add_subset_element!(all_nodes_sub, space_idx, 0, node_idx) # Edges if (i > 1) & (i != cut) # i-1 to i in the npol direction edge_idx1 = edge_start1 + iii * e1_per_i + jj edges[edge_idx1].nodes = [node_idx, node_idx - n_per_i] - SOLPS2IMAS.add_subset_element!(ext_edges_sub, space_idx, 1, edge_idx1) - SOLPS2IMAS.add_subset_element!(ext_xedges_sub, space_idx, 1, edge_idx1) - SOLPS2IMAS.add_subset_element!(all_edges_sub, space_idx, 1, edge_idx1) - SOLPS2IMAS.add_subset_element!(all_xedges_sub, space_idx, 1, edge_idx1) + add_subset_element!(ext_edges_sub, space_idx, 1, edge_idx1) + add_subset_element!(ext_xedges_sub, space_idx, 1, edge_idx1) + add_subset_element!(all_edges_sub, space_idx, 1, edge_idx1) + add_subset_element!(all_xedges_sub, space_idx, 1, edge_idx1) end if (j > 1) # j-1 to j in the nlvl direction edge_idx2 = edge_start2 + ii * e2_per_i + jjj edges[edge_idx2].nodes = [node_idx, node_idx - n_per_j] - SOLPS2IMAS.add_subset_element!(ext_edges_sub, space_idx, 1, edge_idx2) - SOLPS2IMAS.add_subset_element!(ext_yedges_sub, space_idx, 1, edge_idx2) - SOLPS2IMAS.add_subset_element!(all_edges_sub, space_idx, 1, edge_idx2) - SOLPS2IMAS.add_subset_element!(all_yedges_sub, space_idx, 1, edge_idx2) + add_subset_element!(ext_edges_sub, space_idx, 1, edge_idx2) + add_subset_element!(ext_yedges_sub, space_idx, 1, edge_idx2) + add_subset_element!(all_edges_sub, space_idx, 1, edge_idx2) + add_subset_element!(all_yedges_sub, space_idx, 1, edge_idx2) end # Cells @@ -811,8 +812,8 @@ function record_regular_mesh!( node_idx - n_per_i - n_per_j, node_idx - n_per_j, ] - SOLPS2IMAS.add_subset_element!(ext_cells_sub, space_idx, 2, cell_idx) - SOLPS2IMAS.add_subset_element!(all_cells_sub, space_idx, 2, cell_idx) + add_subset_element!(ext_cells_sub, space_idx, 2, cell_idx) + add_subset_element!(all_cells_sub, space_idx, 2, cell_idx) end end end @@ -910,7 +911,7 @@ function cached_mesh_extension!( data["z"] = mesh_z # YAML.write_file(cached_ext_name, data) open(cached_ext_name, "w") do f - JSON.print(f, data) + return JSON.print(f, data) end # fr = open("mesh_r.dat", "w") # fz = open("mesh_z.dat", "w") diff --git a/test/runtests.jl b/test/runtests.jl index 0547680..f75d0a9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,7 @@ using Test using Unitful: Unitful using Interpolations: Interpolations using ArgParse: ArgParse +using GGDUtils: GGDUtils, get_grid_subset_with_index function parse_commandline() s = ArgParse.ArgParseSettings(; description="Run tests. Default is all tests.") @@ -230,7 +231,7 @@ if args["edge_profile_extension"] grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] extended_subs = 1:5 orig_subs = [ - deepcopy(SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, i)) for + deepcopy(get_grid_subset_with_index(grid_ggd, i)) for i ∈ extended_subs ] cfn = SD4SOLPS.cached_mesh_extension!(dd, eqdsk, b2fgmtry; clear_cache=true) @@ -238,9 +239,9 @@ if args["edge_profile_extension"] SD4SOLPS.cached_mesh_extension!(dd, eqdsk, b2fgmtry; grid_ggd_idx=grid_ggd_idx) for j ∈ extended_subs orig_sub = orig_subs[j] - std_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -j) - all_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, j) - ext_sub = SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, -200 - j) + std_sub = get_grid_subset_with_index(grid_ggd, -j) + all_sub = get_grid_subset_with_index(grid_ggd, j) + ext_sub = get_grid_subset_with_index(grid_ggd, -200 - j) no = length(orig_sub.element) ns = length(std_sub.element) na = length(all_sub.element) From 41c51d2da800d9e9a5b51a131003e62978425ed5 Mon Sep 17 00:00:00 2001 From: David Eldon Date: Tue, 7 Nov 2023 12:32:34 -0800 Subject: [PATCH 08/12] Populate boundaries for extended mesh --- src/supersize_profile.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index bb237a3..b1c9235 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -789,6 +789,9 @@ function record_regular_mesh!( if (i > 1) & (i != cut) # i-1 to i in the npol direction edge_idx1 = edge_start1 + iii * e1_per_i + jj edges[edge_idx1].nodes = [node_idx, node_idx - n_per_i] + resize!(edges[edge_idx1].boundary, 2) + edges[edge_idx1].boundary[1].index = node_idx + edges[edge_idx1].boundary[2].index = node_idx - n_per_i add_subset_element!(ext_edges_sub, space_idx, 1, edge_idx1) add_subset_element!(ext_xedges_sub, space_idx, 1, edge_idx1) add_subset_element!(all_edges_sub, space_idx, 1, edge_idx1) @@ -797,6 +800,9 @@ function record_regular_mesh!( if (j > 1) # j-1 to j in the nlvl direction edge_idx2 = edge_start2 + ii * e2_per_i + jjj edges[edge_idx2].nodes = [node_idx, node_idx - n_per_j] + resize!(edges[edge_idx2].boundary, 2) + edges[edge_idx2].boundary[1].index = node_idx + edges[edge_idx2].boundary[2].index = node_idx - n_per_j add_subset_element!(ext_edges_sub, space_idx, 1, edge_idx2) add_subset_element!(ext_yedges_sub, space_idx, 1, edge_idx2) add_subset_element!(all_edges_sub, space_idx, 1, edge_idx2) @@ -812,6 +818,12 @@ function record_regular_mesh!( node_idx - n_per_i - n_per_j, node_idx - n_per_j, ] + resize!(cells[cell_idx].boundary, 4) + cells[cell_idx].boundary[1].index = edge_idx1 + cells[cell_idx].boundary[2].index = edge_idx2 + cells[cell_idx].boundary[3].index = edge_idx1 - 1 + cells[cell_idx].boundary[4].index = edge_idx2 - e2_per_i + add_subset_element!(ext_cells_sub, space_idx, 2, cell_idx) add_subset_element!(all_cells_sub, space_idx, 2, cell_idx) end From 6f838dad37b9aa48c57bdf5ad734490d0b00153d Mon Sep 17 00:00:00 2001 From: David Eldon Date: Tue, 7 Nov 2023 16:09:22 -0800 Subject: [PATCH 09/12] Fixup in response to review - Fix indexing of spaces - Clean up loops - Reduce some function interfaces --- src/supersize_profile.jl | 96 ++++++++++++++++++---------------------- test/runtests.jl | 12 ++--- 2 files changed, 48 insertions(+), 60 deletions(-) diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index b1c9235..cb75739 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -6,9 +6,9 @@ Utilities for extrapolating profiles using OMAS: OMAS using Interpolations: Interpolations using GGDUtils: - GGDUtils, get_grid_subset_with_index, add_subset_element!, get_subset_boundary + GGDUtils, get_grid_subset_with_index, add_subset_element!, get_subset_boundary, + project_prop_on_subset!, get_subset_centers using PolygonOps: PolygonOps -using SOLPS2IMAS: SOLPS2IMAS using JSON: JSON export extrapolate_core @@ -173,7 +173,7 @@ function fill_in_extrapolated_core_profile!( quantity_str = getproperty(quantity_str, Symbol(tag)) end - midplane_cell_centers, quantity = GGDUtils.project_prop_on_subset!( + midplane_cell_centers, quantity = project_prop_on_subset!( quantity_str, cell_subset, midplane_subset, @@ -378,7 +378,7 @@ function mesh_psi_spacing( space = grid_ggd.space[space_idx] midplane_subset = get_grid_subset_with_index(grid_ggd, 11) - midplane_cell_centers = GGDUtils.get_subset_centers(space, midplane_subset) + midplane_cell_centers = get_subset_centers(space, midplane_subset) r_mesh = [midplane_cell_centers[i][1] for i ∈ eachindex(midplane_cell_centers)] z_mesh = [midplane_cell_centers[i][2] for i ∈ eachindex(midplane_cell_centers)] psin_mesh = rzpi.(r_mesh, z_mesh) @@ -429,7 +429,7 @@ function pick_extension_psi_range( grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] space = grid_ggd.space[space_idx] midplane_subset = get_grid_subset_with_index(grid_ggd, 11) - midplane_cell_centers = GGDUtils.get_subset_centers(space, midplane_subset) + midplane_cell_centers = get_subset_centers(space, midplane_subset) psin_midplane = rzpi.(midplane_cell_centers[end][1], midplane_cell_centers[end][2]) # Choose contour levels @@ -469,14 +469,7 @@ grid_ggd_idx: index within ggd space_idx: space number / index of the space to work with within edge_profiles Returns a tuple with vectors of R and Z starting points. """ -function pick_mesh_ext_starting_points( - dd::OMAS.dd; - grid_ggd_idx::Int64=1, - space_idx::Int64=1, -) - grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] - space = grid_ggd.space[space_idx] - +function pick_mesh_ext_starting_points(grid_ggd, space) # Choose starting points for the orthogonal (to the contour) gridlines # Use the existing cells of the standard mesh all_cell_subset = get_grid_subset_with_index(grid_ggd, 5) @@ -484,13 +477,11 @@ function pick_mesh_ext_starting_points( core_edges = get_grid_subset_with_index(grid_ggd, 15) outer_target = get_grid_subset_with_index(grid_ggd, 13) inner_target = get_grid_subset_with_index(grid_ggd, 14) - ci = [core_edges.element[i].object[1].index for i ∈ 1:length(core_edges.element)] - oi = - [outer_target.element[i].object[1].index for i ∈ 1:length(outer_target.element)] - ii = - [inner_target.element[i].object[1].index for i ∈ 1:length(inner_target.element)] + ci = [ele.object[1].index for ele ∈ core_edges.element] + oi = [ele.object[1].index for ele ∈ outer_target.element] + ii = [ele.object[1].index for ele ∈ inner_target.element] border_edges = [] - for i ∈ 1:length(all_border_edges) + for i ∈ eachindex(all_border_edges) bi = all_border_edges[i].object[1].index if !(bi in oi) & !(bi in ii) & !(bi in ci) border_edges = [border_edges; all_border_edges[i]] @@ -522,13 +513,13 @@ of psi_N is large compared to its gradient). r_eq: Equilibrium reconstruction's grid, R coordinates z_eq: Equilibrium reconstruction's grid, Z coordinates psin_eq: Normalized poloidal flux in the equilibrium reconstruction as a function of R and Z -rzpi: linear interpolation of psin_eq() as a function of r_eq and z_eq - This was probably already computed and I think time would be saved by reusing it. - If you don't already have it, you can pass in nothing and let this function calculate it. rstart: R coordinates of starting points for the gradient following. zstart: Z coordinates of starting points nlvl: number of nodes to drop while following the gradient dpsin: node spacing in delta psi_N +rzpi: linear interpolation of psin_eq() as a function of r_eq and z_eq + This was probably already computed and I think time would be saved by reusing it. + If you don't already have it, you can pass in nothing and let this function calculate it. Returns two matrices with R and Z coordinates of the mesh extension """ @@ -537,11 +528,11 @@ function mesh_ext_follow_grad( r_eq::Vector{Float64}, z_eq::Vector{Float64}, psin_eq::Matrix, - rzpi, rstart::Vector{Float64}, zstart::Vector{Float64}, nlvl::Int64, dpsin::Float64, + rzpi=nothing, ) npol = length(rstart) mesh_r = zeros((npol, nlvl)) @@ -598,7 +589,11 @@ eqt: equilibrium.time_slice information mesh_r: matrix of R values for the extended mesh mesh_z: matrix of Z values for the extended mesh """ -function modify_mesh_ext_near_x!(eqt, mesh_r, mesh_z) +function modify_mesh_ext_near_x!( + eqt::OMAS.equilibrium__time_slice, + mesh_r::Matrix{Float64}, + mesh_z::Matrix{Float64}, +) # There's a special path; the one that probably should've gone through the # secondary X-point (if there is one). Any numerical error will make this # path miss the X-point and go off to somewhere crazy instead, so instead of @@ -681,9 +676,8 @@ end record_regular_mesh!() Records arrays of mesh data from regular 2D arrays into the DD -dd: the data dictionary -grid_ggd_idx: index of the grid_ggd within edge_profiles -space_idx: index of the space in edge_profiles +grid_ggd: grid_ggd within edge_profiles +space: space in edge_profiles mesh_r: Matrix of R values along a mesh. Should be 2D. The two dimensions are in the radial and poloidal directions. mesh_z: Z values to go with mesh_r. @@ -693,20 +687,19 @@ cut: Poloidal index of a cut between two groups of disconnected cells. Poloidal """ #!format on function record_regular_mesh!( - dd::OMAS.dd, - grid_ggd_idx::Int64, - space_idx::Int64, + grid_ggd, + space, mesh_r::Matrix{Float64}, mesh_z::Matrix{Float64}, cut::Int64, ) - grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] - space = grid_ggd.space[space_idx] - npol, nlvl = size(mesh_r) o0 = space.objects_per_dimension[1] # Nodes o1 = space.objects_per_dimension[2] # Edges o2 = space.objects_per_dimension[3] # Cells (2D projection) + node_dim = 1 + edge_dim = 2 + cell_dim = 3 n_old_node = length(o0.object) n_new_node = nlvl * npol @@ -770,6 +763,7 @@ function record_regular_mesh!( nodes = resize!(o0.object, n_nodes) edges = resize!(o1.object, n_edges) cells = resize!(o2.object, n_cells) + space_idx = space.identifier.index for i ∈ 1:npol pastpc = i > cut for j ∈ 1:nlvl @@ -782,8 +776,8 @@ function record_regular_mesh!( # Nodes node_idx = node_start + ii * n_per_i + jj nodes[node_idx].geometry = [mesh_r[i, j], mesh_z[i, j]] - add_subset_element!(ext_nodes_sub, space_idx, 0, node_idx) - add_subset_element!(all_nodes_sub, space_idx, 0, node_idx) + add_subset_element!(ext_nodes_sub, space_idx, node_dim, node_idx) + add_subset_element!(all_nodes_sub, space_idx, node_dim, node_idx) # Edges if (i > 1) & (i != cut) # i-1 to i in the npol direction @@ -792,10 +786,10 @@ function record_regular_mesh!( resize!(edges[edge_idx1].boundary, 2) edges[edge_idx1].boundary[1].index = node_idx edges[edge_idx1].boundary[2].index = node_idx - n_per_i - add_subset_element!(ext_edges_sub, space_idx, 1, edge_idx1) - add_subset_element!(ext_xedges_sub, space_idx, 1, edge_idx1) - add_subset_element!(all_edges_sub, space_idx, 1, edge_idx1) - add_subset_element!(all_xedges_sub, space_idx, 1, edge_idx1) + add_subset_element!(ext_edges_sub, space_idx, edge_dim, edge_idx1) + add_subset_element!(ext_xedges_sub, space_idx, edge_dim, edge_idx1) + add_subset_element!(all_edges_sub, space_idx, edge_dim, edge_idx1) + add_subset_element!(all_xedges_sub, space_idx, edge_dim, edge_idx1) end if (j > 1) # j-1 to j in the nlvl direction edge_idx2 = edge_start2 + ii * e2_per_i + jjj @@ -803,10 +797,10 @@ function record_regular_mesh!( resize!(edges[edge_idx2].boundary, 2) edges[edge_idx2].boundary[1].index = node_idx edges[edge_idx2].boundary[2].index = node_idx - n_per_j - add_subset_element!(ext_edges_sub, space_idx, 1, edge_idx2) - add_subset_element!(ext_yedges_sub, space_idx, 1, edge_idx2) - add_subset_element!(all_edges_sub, space_idx, 1, edge_idx2) - add_subset_element!(all_yedges_sub, space_idx, 1, edge_idx2) + add_subset_element!(ext_edges_sub, space_idx, edge_dim, edge_idx2) + add_subset_element!(ext_yedges_sub, space_idx, edge_dim, edge_idx2) + add_subset_element!(all_edges_sub, space_idx, edge_dim, edge_idx2) + add_subset_element!(all_yedges_sub, space_idx, edge_dim, edge_idx2) end # Cells @@ -824,8 +818,8 @@ function record_regular_mesh!( cells[cell_idx].boundary[3].index = edge_idx1 - 1 cells[cell_idx].boundary[4].index = edge_idx2 - e2_per_i - add_subset_element!(ext_cells_sub, space_idx, 2, cell_idx) - add_subset_element!(all_cells_sub, space_idx, 2, cell_idx) + add_subset_element!(ext_cells_sub, space_idx, cell_dim, cell_idx) + add_subset_element!(all_cells_sub, space_idx, cell_dim, cell_idx) end end end @@ -899,9 +893,8 @@ function cached_mesh_extension!( npol = md["npol"] nlvl = md["nlvl"] record_regular_mesh!( - dd, - grid_ggd_idx, - space_idx, + dd.edge_profiles.grid_ggd[grid_ggd_idx], + dd.edge_profiles.grid_ggd[grid_ggd_idx].space[space_idx], mesh_r, mesh_z, pfr_transition, @@ -965,18 +958,17 @@ function mesh_extension_sol!( ) nlvl = length(psin_levels) dpsin = psin_levels[2] - psin_levels[1] - grad_start_r, grad_start_z = - pick_mesh_ext_starting_points(dd; grid_ggd_idx, space_idx) + grad_start_r, grad_start_z = pick_mesh_ext_starting_points(grid_ggd, space) npol = length(grad_start_r) mesh_r, mesh_z = mesh_ext_follow_grad( r_eq, z_eq, psin_eq, - rzpi, grad_start_r, grad_start_z, nlvl, dpsin, + rzpi, ) modify_mesh_ext_near_x!(eqt, mesh_r, mesh_z) @@ -987,7 +979,7 @@ function mesh_extension_sol!( # The PFR flag needs to be respected; pfr nodes don't connect to non-pfr nodes pfr = rzpi.(mesh_r[:, 1], mesh_z[:, 1]) .< 1 pfr_transition = argmax(abs.(diff(pfr))) - record_regular_mesh!(dd, grid_ggd_idx, space_idx, mesh_r, mesh_z, pfr_transition) + record_regular_mesh!(grid_ggd, space, mesh_r, mesh_z, pfr_transition) return mesh_r, mesh_z, pfr_transition end diff --git a/test/runtests.jl b/test/runtests.jl index f75d0a9..ff6f825 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -242,14 +242,10 @@ if args["edge_profile_extension"] std_sub = get_grid_subset_with_index(grid_ggd, -j) all_sub = get_grid_subset_with_index(grid_ggd, j) ext_sub = get_grid_subset_with_index(grid_ggd, -200 - j) - no = length(orig_sub.element) - ns = length(std_sub.element) - na = length(all_sub.element) - ne = length(ext_sub.element) - orig_indices = [orig_sub.element[i].object[1].index for i ∈ 1:no] - std_indices = [std_sub.element[i].object[1].index for i ∈ 1:ns] - all_indices = [all_sub.element[i].object[1].index for i ∈ 1:na] - ext_indices = [ext_sub.element[i].object[1].index for i ∈ 1:ne] + orig_indices = [ele.object[1].index for ele ∈ orig_sub.element] + std_indices = [ele.object[1].index for ele ∈ std_sub.element] + all_indices = [ele.object[1].index for ele ∈ all_sub.element] + ext_indices = [ele.object[1].index for ele ∈ ext_sub.element] @test std_sub.identifier.index == -j @test all(orig_indices .== std_indices) all_indices_reconstruct = [std_indices; ext_indices] From b81facb1ae0a3b7e1289983c46d3587c93a06689 Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Wed, 17 Jan 2024 11:51:03 -0800 Subject: [PATCH 10/12] Rolling out v.0.2.0 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f727ad2..a3e1890 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SD4SOLPS" uuid = "d8db6f1b-e564-4c04-bba3-ef399f78c070" authors = ["David Eldon "] -version = "0.1.0" +version = "0.2.0" [deps] ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" From 120e2406478fff3e3d42021962557c5329b7dfb8 Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Wed, 17 Jan 2024 18:43:50 -0800 Subject: [PATCH 11/12] Adding a demonstration of whole project --- example/demo.ipynb | 461 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 461 insertions(+) create mode 100644 example/demo.ipynb diff --git a/example/demo.ipynb b/example/demo.ipynb new file mode 100644 index 0000000..6841ccf --- /dev/null +++ b/example/demo.ipynb @@ -0,0 +1,461 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using Pkg;\n", + "Pkg.activate(\"./\")\n", + "Pkg.instantiate()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Loading SOLPS output into IMAS data structure" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using SOLPS2IMAS: SOLPS2IMAS" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "b2gmtry = \"samples/b2fgmtry\"\n", + "b2output = \"samples/b2time.nc\"\n", + "gsdesc = \"samples/gridspacedesc.yml\"\n", + "b2mn = \"samples/b2mn.dat\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ids = SOLPS2IMAS.solps2imas(b2gmtry, b2output, gsdesc, b2mn)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Visualising some properties" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using GGDUtils: GGDUtils\n", + "using Plots" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "IMAS ids store mesh information for edge profiles in grid_ggd. There are options to have multiple space representaions but typically you will have only one space describing the SOLPS mesh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "grid_ggd = ids.edge_profiles.grid_ggd[1] # First grid_ggd time slice. It is allowed to vary in time\n", + "space = grid_ggd.space[1] # First space in this grid_ggd" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plotting grid and subsets" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(space) # Simply plot the grid described in space, all common arguments to plot can be given here\n", + "\n", + "# You can overlay any subset by giving a second argument\n", + "# Labels \n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 6), markercolor=:chocolate1)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 7), linecolor=:red, linewidth=2)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 8), linecolor=:darkred, linewidth=2)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 9), linecolor=:limegreen, linewidth=2)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 10), linecolor=:darkgreen, linewidth=2)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 11), linecolor=:cyan, linewidth=2)\n", + "# plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 12), linecolor=:teal, linewidth=1)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 13), linecolor=:royalblue1, linewidth=2)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 14), linecolor=:navyblue, linewidth=2)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 15), linecolor=:fuchsia, linewidth=2, linestyle=:dash)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 16), linecolor=:purple4, linewidth=2, linestyle=:dash)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 101), markershape=:rect, markercolor=:royalblue1)\n", + "# plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 102), markershape=:rect, markercolor=:maroon)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 103), markershape=:diamond, markercolor=:fuchsia)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 104), markershape=:diamond, markercolor=:purple4)\n", + "\n", + "# Legend is supressed unless asked for specifically\n", + "plot!(legend=true)\n", + "# Default labels are subset.identifier.name but can be changed by providing a label argument\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plotting 2D quantities as heatmaps" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "n_e = GGDUtils.get_prop_with_grid_subset_index(ids.edge_profiles.ggd[1].electrons.density, 5)\n", + "plot(ids.edge_profiles.grid_ggd, n_e, colorbar_title=\"Electrons density / m^(-3)\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can overlap any grid on top of a quantity" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(ids.edge_profiles.grid_ggd, n_e) # Note default label in colorbar\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 16), linecolor=:black, linewidth=2, linestyle=:solid, label=\"Separatix\", legend=true)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Adding equilibrium data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import SD4SOLPS: SD4SOLPS" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "eqdsk = \"samples/g002296.00200\"\n", + "SD4SOLPS.geqdsk_to_imas!(eqdsk, ids)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Extrapole core profile" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "SD4SOLPS.fill_in_extrapolated_core_profile!(ids, \"electrons.density\"; method=\"simple\")\n", + "SD4SOLPS.fill_in_extrapolated_core_profile!(ids, \"electrons.temperature\"; method=\"simple\")\n", + "# ... more profiles here as they become available in b2time" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Extend mesh outside SOLPS mesh to the device wall" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "SD4SOLPS.cached_mesh_extension!(ids, eqdsk, b2gmtry)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Loading a synthetic diagnostic" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using SynthDiag: SynthDiag" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Add interferometer chord details using a json file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "SynthDiag.add_interferometer!(\"samples/default_interferometer.json\", ids)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plotting the interferometer geometry on top of SOLPS mesh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(space)\n", + "plot!(ids.interferometer) # Default plot_type is :los \n", + "plot!(legend=true)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can provide custom length and thickness of mirror to be plotted and linewidth of the laser beams" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(space)\n", + "plot!(ids.interferometer, mirror_length=0.7, linewidth=4, mirror_thickness=0.2)\n", + "plot!(legend=true)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Or you can choose to omit the mirror" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(space)\n", + "plot!(ids.interferometer, mirror=false)\n", + "plot!(legend=true)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can plot a single channel as well. You can override the in-built channel name for the label." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(space)\n", + "plot!(ids.interferometer.channel[1], label=\"Channel 1\")\n", + "plot!(legend=true)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plotting interferometer data vs time\n", + "\n", + " * Use plot_type=:n_e for integrated electron density data\n", + " * Use plot_type=:n_e_average for averaged electron density data\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(ids.interferometer, plot_type=:n_e)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(ids.interferometer, plot_type=:n_e_average)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Again, to plot an individual channel, just provide the channel with correct plot_type" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(ids.interferometer.channel[1], plot_type=:n_e_average)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Load Langmuir Probes\n", + "\n", + "Same syntax as interferometer" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "SynthDiag.add_langmuir_probes!(\"samples/default_langmuir_probes.json\", ids)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " Data visualization recipes for langmuir probes have not been created yet and might not be created but you can still see them using built in plotting methodsan d knowledgeof IMAS data structure." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(ids.langmuir_probes.embedded[1].time, ids.langmuir_probes.embedded[1].n_e.data, label=ids.langmuir_probes.embedded[1].name)\n", + "plot!(ylabel=\"Electron density / m^(-3)\", xlabel=\"Time / s\", legend=true)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.9.2", + "language": "julia", + "name": "julia-1.9" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.9.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From 48d6386e068dcea185abdd14bb4b92a7ffb1fc4f Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Thu, 18 Jan 2024 16:33:34 -0800 Subject: [PATCH 12/12] Adding makefile and updated README * makefile has been added to allow easier creation of working environment * README has been updated to guide the user through the process of creating a working environment and added links to setup and demonstration wiki. --- README.md | 52 ++++++++++++++++++++++++++++++++ example/create_env_with_clone.sh | 38 +++++++++++++++++++++++ example/create_env_with_git.sh | 11 +++++++ example/makefile | 25 +++++++++++++++ example/rm_cloned_repo.sh | 32 ++++++++++++++++++++ 5 files changed, 158 insertions(+) create mode 100644 example/create_env_with_clone.sh create mode 100644 example/create_env_with_git.sh create mode 100644 example/makefile create mode 100644 example/rm_cloned_repo.sh diff --git a/README.md b/README.md index b4143c7..ed84d70 100644 --- a/README.md +++ b/README.md @@ -9,3 +9,55 @@ Steps: 2) Load equilibrium (that the SOLPS mesh was based on) into IMAS DD format 3) Make assumptions to extend profiles into the core and far SOL, if needed 4) Run synthetic diagnostic models and record output + + +## Installation + +### Using make file + +#### Option 1: Download the [example directory](https://github.com/ProjectTorreyPines/SD4SOLPS.jl/tree/master/example) in this repo: + +```bash +% make help +Help Menu + +make env_with_cloned_repo: Creates a Julia environment with the cloned repositories +make env_with_git_url: Creates a Julia environment with the git urls without creating local clones +make clean: Deletes Project.toml and Manifest.toml for a fresh start +make Clean: Deletes Project.toml, Manifest.toml, and any cloned repositories for a fresh start +``` +##### Environment with cloned repos + +`make env_with_cloned_repo` will clone all required git repos one level up from your current working directory and use them to define the environment files that will be generated in your working directory. This way, in future you can do `git pull` or change branches to update your working environment. + +##### Static environment with git url + +`make env_with_git_url` will not create local git clones for you, but it will still clone master branches of all the required repos in julia package management directory and will create a static environment for you that will remain tied to the version of packages when this environment was created. + +##### Clean up + +`make clean` or `make Clean` can be used clean up environment files and redo the environment generation. + +#### Option 2: Clone this repo + +```bash +% git clone git@github.com:ProjectTorreyPines/SD4SOLPS.jl.git +% cd example +$ make help +Help Menu + +make env_with_cloned_repo: Creates a Julia environment with the cloned repositories +make env_with_git_url: Creates a Julia environment with the git urls without creating local clones +make clean: Deletes Project.toml and Manifest.toml for a fresh start +make Clean: Deletes Project.toml, Manifest.toml, and any cloned repositories for a fresh start +``` + +Further options are same as above except for the difference that in case of cloning local copies of repos, they will be kept on same level as where you cloned SD4SOLPS.jl repo. + +### Manual Install + +Refer to the instructions on this [wiki page](https://github.com/ProjectTorreyPines/SD4SOLPS.jl/wiki). + +## Use + +Refer to the instructions on this [wiki page](https://github.com/ProjectTorreyPines/SD4SOLPS.jl/wiki/Demo). diff --git a/example/create_env_with_clone.sh b/example/create_env_with_clone.sh new file mode 100644 index 0000000..f1c4225 --- /dev/null +++ b/example/create_env_with_clone.sh @@ -0,0 +1,38 @@ +#!/bin/bash + +is_repo=$(git rev-parse --is-inside-work-tree 2>/dev/null) +if [[ ${is_repo} == "true" ]] +then + SD4SOLPS_path="$(git rev-parse --show-toplevel)" + project_path="$(dirname ${SD4SOLPS_path})" +else + SD4SOLPS_path="" + project_path="$(dirname $(pwd))" +fi +echo "Cloning git repositories to "$project_path +orig_path=$(pwd) +cd $project_path +if [[ -z $SD4SOLPS_path ]] +then + git clone git@github.com:ProjectTorreyPines/SD4SOLPS.jl.git + SD4SOLPS_path=$project_path"/SD4SOLPS.jl" +fi +# in a for loop clone SynthDiag.jl SOLPS2IMAS.jl OMAS.jl GGDUtils.jl +for repo in SynthDiag.jl SOLPS2IMAS.jl OMAS.jl GGDUtils.jl +do + if [[ ! -d $repo ]] + then + git clone "git@github.com:ProjectTorreyPines/"$repo".git" + fi +done +cd $orig_path + +echo "Generating Project.toml and Manifest.toml" +OMAS_path=$project_path"/OMAS.jl" +EFIT_url="git@github.com:JuliaFusion/EFIT.jl.git" +julia --project=. -e 'using Pkg; Pkg.develop(path="'$OMAS_path'"); Pkg.instantiate()' +julia --project=. -e 'using Pkg; Pkg.add(url="'$EFIT_url'", rev="master"); Pkg.instantiate()' +for repo in GGDUtils SOLPS2IMAS SynthDiag SD4SOLPS; do + repo_path=$project_path"/"$repo".jl" + julia --project=. -e 'using Pkg; Pkg.develop(path="'${repo_path}'"); Pkg.instantiate()' +done \ No newline at end of file diff --git a/example/create_env_with_git.sh b/example/create_env_with_git.sh new file mode 100644 index 0000000..38188f9 --- /dev/null +++ b/example/create_env_with_git.sh @@ -0,0 +1,11 @@ +#!/bin/bash + +echo "Generating Project.toml and Manifest.toml" +OMAS_url="git@github.com:ProjectTorreyPines/OMAS.jl.git" +EFIT_url="git@github.com:JuliaFusion/EFIT.jl.git" +julia --project=. -e 'using Pkg; Pkg.add(url="'$OMAS_url'", rev="master"); Pkg.instantiate()' +julia --project=. -e 'using Pkg; Pkg.add(url="'$EFIT_url'", rev="master"); Pkg.instantiate()' +for repo in GGDUtils SOLPS2IMAS SynthDiag SD4SOLPS; do + repo_url="git@github.com:ProjectTorreyPines/"$repo".jl.git" + julia --project=. -e 'using Pkg; Pkg.add(url="'$repo_url'", rev="master"); Pkg.instantiate()' +done \ No newline at end of file diff --git a/example/makefile b/example/makefile new file mode 100644 index 0000000..8d6c7a2 --- /dev/null +++ b/example/makefile @@ -0,0 +1,25 @@ +SHELL := /bin/zsh +help: + @echo "Help Menu" + @echo + @echo "make env_with_cloned_repo: Creates a Julia environment with the cloned repositories" + @echo "make env_with_git_url: Creates a Julia environment with the git urls without creating local clones" + @echo "make clean: Deletes Project.toml and Manifest.toml for a fresh start" + @echo "make Clean: Deletes Project.toml, Manifest.toml, and any cloned repositories for a fresh start" + @echo + +env_with_cloned_repo: + chmod +x create_env_with_clone.sh + ./create_env_with_clone.sh + +env_with_git_url: + chmod +x create_env_with_git.sh + ./create_env_with_git.sh + +clean: + @echo "Deleting Project.toml and Manifest.toml" + rm Project.toml Manifest.toml + +Clean: clean + chmod +x rm_cloned_repo.sh + ./rm_cloned_repo.sh \ No newline at end of file diff --git a/example/rm_cloned_repo.sh b/example/rm_cloned_repo.sh new file mode 100644 index 0000000..351e0cc --- /dev/null +++ b/example/rm_cloned_repo.sh @@ -0,0 +1,32 @@ +#!/bin/bash + +echo "Deleting any cloned repo" +is_repo=$(git rev-parse --is-inside-work-tree 2>/dev/null) +if [[ ${is_repo} == "true" ]] +then + SD4SOLPS_path="$(git rev-parse --show-toplevel)" + project_path="$(dirname ${SD4SOLPS_path})" +else + SD4SOLPS_path="" + project_path="$(dirname $(pwd))" +fi + +echo "Removing extra cloned git repositories from "$project_path + +if [[ -z $SD4SOLPS_path ]] +then + # Check if $project_path"/SD4SOLPS.jl" exists + if [[ -d $project_path"/SD4SOLPS.jl" ]] + then + echo "Removing "$project_path"/SD4SOLPS.jl" + rm -rf $project_path"/SD4SOLPS.jl" + fi +fi +for repo in SynthDiag.jl SOLPS2IMAS.jl OMAS.jl GGDUtils.jl +do + if [[ -d $project_path"/"$repo ]] + then + echo "Removing "$project_path"/"$repo + rm -rf $project_path"/"$repo + fi +done \ No newline at end of file