From 58c6c578ee67b561493997aa9c0071fe1c3f519b Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 8 Feb 2023 10:50:27 -0800 Subject: [PATCH 1/4] Improve grid interpolation for restart files - Linear interpolation option - Floor density and pressure to non-zero values --- boutdata/restart.py | 29 ++++++++++++++++++++++++++--- 1 file changed, 26 insertions(+), 3 deletions(-) diff --git a/boutdata/restart.py b/boutdata/restart.py index c4571f8..6074386 100644 --- a/boutdata/restart.py +++ b/boutdata/restart.py @@ -974,6 +974,7 @@ def change_grid( output=".", interpolator="nearest", show=False, + floor=True, ): """ Convert a set of restart files from one grid to another @@ -992,9 +993,11 @@ def change_grid( output : str, optional Directory where output restart files will be written interpolator : str, optional - Interpolation method to use. Options are 'nearest', 'CloughTocher', 'RBF' + Interpolation method to use. Options are 'nearest', 'CloughTocher', 'RBF', 'linear' show : bool, optional Display the interpolated fields using Matplotlib + floor : bool, optional + Apply floor to interpolated densities and pressures """ @@ -1047,7 +1050,10 @@ def change_grid( # Read information from a restart file with DataFile(file_list[0]) as f: for var in copy_vars: - copy_data[var] = f[var] + try: + copy_data[var] = f[var] + except: + print("Failed to copy variable '{}'".format(var)) # Get a list of variables varnames = f.list() @@ -1109,7 +1115,6 @@ def extrapolate_yguards(data2d, myg): from_data.flatten(), neighbors=50, ) - elif interpolator == "nearest": # Nearest neighbour. Tends to be robust from scipy.interpolate import NearestNDInterpolator @@ -1117,6 +1122,13 @@ def extrapolate_yguards(data2d, myg): interp = NearestNDInterpolator( list(zip(from_Rxy.flatten(), from_Zxy.flatten())), from_data.flatten() ) + elif interpolator == "linear": + # Linear interpolator + from scipy.interpolate import LinearNDInterpolator + interp = LinearNDInterpolator( + list(zip(from_Rxy.flatten(), from_Zxy.flatten())), from_data.flatten(), + fill_value = 0.0 + ) else: raise ValueError("Invalid interpolator") @@ -1132,6 +1144,17 @@ def extrapolate_yguards(data2d, myg): np.amax(to_data), ) ) + + if floor: + if var[0] == 'N' and var[1] != 'V': + # A density + to_data = np.clip(to_data, 1e-5, None) + print("\t\t floor -> {}:{}".format(np.amin(to_data), np.amax(to_data))) + elif var[0] == 'P': + # Pressure + to_data = np.clip(to_data, 1e-8, None) + print("\t\t floor -> {}:{}".format(np.amin(to_data), np.amax(to_data))) + if show: import matplotlib.pyplot as plt From 888eb5e7d4eb093d2e321f0272efefaf6178c5b9 Mon Sep 17 00:00:00 2001 From: bendudson Date: Wed, 8 Feb 2023 18:57:00 +0000 Subject: [PATCH 2/4] Apply black changes --- boutdata/mayavi2.py | 1 - boutdata/mms.py | 1 - boutdata/restart.py | 12 ++++++------ 3 files changed, 6 insertions(+), 8 deletions(-) diff --git a/boutdata/mayavi2.py b/boutdata/mayavi2.py index f583388..76a62b3 100644 --- a/boutdata/mayavi2.py +++ b/boutdata/mayavi2.py @@ -44,7 +44,6 @@ def aligned_points(grid, nz=1, period=1.0, maxshift=0.4): def create_grid(grid, data, period=1): - s = np.shape(data) nx = grid["nx"] diff --git a/boutdata/mms.py b/boutdata/mms.py index 7354326..c4d30cf 100644 --- a/boutdata/mms.py +++ b/boutdata/mms.py @@ -408,7 +408,6 @@ def write(self, nx, ny, output, MXG=2): ("sinty", self.sinty), ("zShift", self.zShift), ]: - # Note: This is slow, and could be improved using something like lambdify values = zeros([ngx, ngy]) for i, x in enumerate(xarr): diff --git a/boutdata/restart.py b/boutdata/restart.py index 6074386..5a74dc5 100644 --- a/boutdata/restart.py +++ b/boutdata/restart.py @@ -174,7 +174,6 @@ def is_pow2(x): # Open the restart file in read mode and create the new file with DataFile(f) as old, DataFile(new_f, write=True, create=True) as new: - # Find the dimension for var in old.list(): # Read the data @@ -223,7 +222,6 @@ def is_pow2(x): # Find 3D variables if old.ndims(var) == 3: - # Asynchronous call (locks first at .get()) jobs.append( pool.apply_async( @@ -1125,9 +1123,11 @@ def extrapolate_yguards(data2d, myg): elif interpolator == "linear": # Linear interpolator from scipy.interpolate import LinearNDInterpolator + interp = LinearNDInterpolator( - list(zip(from_Rxy.flatten(), from_Zxy.flatten())), from_data.flatten(), - fill_value = 0.0 + list(zip(from_Rxy.flatten(), from_Zxy.flatten())), + from_data.flatten(), + fill_value=0.0, ) else: raise ValueError("Invalid interpolator") @@ -1146,11 +1146,11 @@ def extrapolate_yguards(data2d, myg): ) if floor: - if var[0] == 'N' and var[1] != 'V': + if var[0] == "N" and var[1] != "V": # A density to_data = np.clip(to_data, 1e-5, None) print("\t\t floor -> {}:{}".format(np.amin(to_data), np.amax(to_data))) - elif var[0] == 'P': + elif var[0] == "P": # Pressure to_data = np.clip(to_data, 1e-8, None) print("\t\t floor -> {}:{}".format(np.amin(to_data), np.amax(to_data))) From 7d20da418554dc87ae014402d4ea999ec33445de Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 5 Jul 2024 09:57:40 -0700 Subject: [PATCH 3/4] Update boutdata/restart.py Co-authored-by: David Bold --- boutdata/restart.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/boutdata/restart.py b/boutdata/restart.py index 5a74dc5..f4a6e4f 100644 --- a/boutdata/restart.py +++ b/boutdata/restart.py @@ -1048,10 +1048,10 @@ def change_grid( # Read information from a restart file with DataFile(file_list[0]) as f: for var in copy_vars: - try: - copy_data[var] = f[var] - except: - print("Failed to copy variable '{}'".format(var)) + if var not in f: + print(f"{var} not present - skipping") + continue + copy_data[var] = f[var] # Get a list of variables varnames = f.list() From 96d1ad6bb6d9c463d6b39753b96e7509b3ba1e0d Mon Sep 17 00:00:00 2001 From: bendudson Date: Fri, 5 Jul 2024 16:58:04 +0000 Subject: [PATCH 4/4] Apply black changes --- boutdata/data.py | 6 +++--- boutdata/gen_surface.py | 1 + 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/boutdata/data.py b/boutdata/data.py index 25ca319..a6a69fb 100644 --- a/boutdata/data.py +++ b/boutdata/data.py @@ -695,9 +695,9 @@ def __init__( comments = [] if inline_comment is not None: parent_section.inline_comments[sectionname] = inline_comment - parent_section._comment_whitespace[ - sectionname - ] = comment_whitespace + parent_section._comment_whitespace[sectionname] = ( + comment_whitespace + ) else: # A key=value pair diff --git a/boutdata/gen_surface.py b/boutdata/gen_surface.py index ed2ce26..808ced9 100644 --- a/boutdata/gen_surface.py +++ b/boutdata/gen_surface.py @@ -1,6 +1,7 @@ """Flux surface generator for tokamak grid files """ + from __future__ import print_function import numpy as np