Added tutorial cases

Examples/Artemis_Tutorial_Cases/README.md

@@ -0,0 +1,18 @@
+# Test Cases
+This file contains two test cases: an air-filled waveguide and a magnetically tunable filter. The simulation parameters are defined in an input file, "/inputs_*". To simulate the air-filled waveguide using NERSC's Perlmutter, run
+> sbatch perl.run
+The simulation will generate a folder of diagnostic data titled "/diags". Next, to convert the diagnostic data into readable txt files, run
+> python3 batch_convert.py 
+The converted data will be stored in "/raw_data*". Finally, to extract the phase constant and attenuation spectrums, run 
+> python3 alphaExtraction.py
+
+This folders contains the following files:
+- alphaExtraction.py : Extracts the components of the wavevector, k, from a "/raw_data*" file
+- batch_convert.py : Converts diagnostic data, "/diags", into readable data, "raw_data*".
+- inputs_3d_empty_X_band : Input file defining simulation parameters for air-filled waveguide case.
+- inputs_3d_LLG_filter : Input file defining simulation parameters for the magnetically tunable filter.
+- main3d.gnu.TPROF.MTMPI.CUDA.ex : Executable file used for air-filled waveguide; the non-LLG case.
+- main3d.gnu.TPROF.MTMPI.CUDA.LLG.ex : Executable file used for the magnetically tunable filter; the LLG case.
+- perl_empty.run : Run file for NERSC's Perlmutter for the non-LLG case.
+- perl_llg.run : Run file for NERSC's Perlmutter for LLG case.
+

Examples/Artemis_Tutorial_Cases/alphaExtraction.py

@@ -0,0 +1,184 @@
+import glob
+import math
+import re
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.interpolate import interp1d
+from scipy.linalg import hankel, svd
+
+
+def esprit(x, p, n_sources):
+    """1D Total Least Squares ESPRIT (TLS-ESPRIT) algorithm."""
+    N = len(x)
+    H = hankel(x[: N - p + 1], x[N - p :])
+
+    U, _, _ = svd(H, full_matrices=False)
+    Us = U[:, :n_sources]
+
+    Us1, Us2 = Us[:-1, :], Us[1:, :]
+
+    K = np.hstack((Us1, Us2))
+    _, _, Vh_k = svd(K)
+    V = Vh_k.conj().T
+
+    V12 = V[:n_sources, n_sources:]
+    V22 = V[n_sources:, n_sources:]
+    Psi = -V12 @ np.linalg.inv(V22)
+
+    z = np.linalg.eigvals(Psi)
+    wm = np.log(z)
+
+    return wm
+
+
+# --- Parameters & Grid Setup ---
+# Physical Constants & Waveguide Setup
+EP0 = 8.854e-12
+MU0 = 4 * math.pi * 1e-7
+CC = 3e8
+
+# GENERAL PARAMETERS
+freq = 10.5e9  # frequency
+XX, YY, ZZ = 14.95e-3, 10.16e-3, 500.0e-3  # dimensions
+NX, NY, NZ = 128, 4, 128  # number of cells
+
+DX, DY, DZ = XX / NX, YY / NY, ZZ / NZ
+
+DT = 3.11275047e-13  # time step in seconds
+TP = 1.0 / freq
+Tstep_end = 1000000  # maximum time steps
+Tstep_int = 100  # interval between time steps
+
+file_pattern = "./raw_data_empty/*.txt"  # data location
+
+k0 = 1
+k1 = NZ - 1
+dk = 1
+NZ_record = int(math.floor(k1 / dk))
+
+Tstep_record = math.floor((4.7 / 6.3 * Tstep_end) / Tstep_int)
+Tstep_trunc = Tstep_record
+
+Time_axis = np.arange(Tstep_int, Tstep_end + 1, Tstep_int)
+
+# Spectral Setup
+fs = 1.0 / (Tstep_int * DT)
+Tstep_zeros = 0
+Tstep_ext = Tstep_record + Tstep_zeros
+NFFT_ext = 1 << (Tstep_ext - 1).bit_length()
+
+f_ext = (fs / 2.0) * np.linspace(0, 1, (NFFT_ext // 2) + 1)
+f_ext_GHz = f_ext / 1e9
+f_ext_step = np.mean(np.diff(f_ext))
+
+# --- Data Import ---
+txt_files = sorted(
+    glob.glob(file_pattern),
+    key=lambda x: int(re.search(r"(\d+)\.txt$", x).group(1)) if re.search(r"(\d+)\.txt$", x) else 0,
+)
+
+Ey_obs_ret = np.zeros((NZ_record, Tstep_record))
+for tcount, fname in enumerate(txt_files[:Tstep_record]):
+    data = np.loadtxt(fname, skiprows=0)
+    Ey_obs_ret[:, tcount] = data[:NZ_record, 1]
+
+ZZ_axis = DZ * np.arange(k0, k1 + 1, dk)
+
+# --- FFT Processing ---
+ey_ext_rec = np.hstack(
+    (
+        Ey_obs_ret[:, :Tstep_trunc],
+        np.zeros((NZ_record, Tstep_record - Tstep_trunc)),
+        np.zeros((NZ_record, Tstep_zeros)),
+    )
+)
+
+ey_extff_rec = np.fft.fft(ey_ext_rec, n=NFFT_ext, axis=1) / Tstep_ext
+
+# --- ESPRIT Processing (Integrated into Main Script) ---
+ko1_m = int(math.floor((ZZ - DZ) / (DZ * dk)))
+ko2_m = int(math.floor(DZ / (DZ * dk)))
+
+num_freqs = int(NFFT_ext // 2) + 1
+beta_esprit = np.zeros(num_freqs)
+alpha_esprit = np.zeros(num_freqs)
+spatial_step = DZ * dk
+
+for fpoint in range(num_freqs):
+    slice_m = ey_extff_rec[ko2_m:ko1_m, fpoint]
+    wm = esprit(slice_m, p=3, n_sources=1)
+
+    beta_esprit[fpoint] = np.abs(np.imag(wm[0])) / spatial_step
+    alpha_esprit[fpoint] = np.abs(np.real(wm[0])) / spatial_step
+
+N_m = len(beta_esprit)
+
+# --- Data Exports ---
+alpha_m_db = 10 * np.log10(np.exp(2)) * (alpha_esprit / 100)
+
+alpha_data = np.column_stack((f_ext_GHz[:N_m], alpha_m_db))
+np.savetxt("alpha.txt", alpha_data, header="Frequency Forward", fmt="%.6f", comments="")
+
+beta_data = np.column_stack((f_ext_GHz[:N_m], beta_esprit))
+np.savetxt("beta.txt", beta_data, header="Frequency Forward", fmt="%.6f", comments="")
+
+Nf_ey = int(round(10.2e9 / f_ext_step))
+ey_data = np.column_stack((ZZ_axis, np.abs(ey_extff_rec[:, Nf_ey])))
+np.savetxt("Ey.txt", ey_data, header="Z position |Ey|", fmt="%.6f", comments="")
+
+# --- Target Plots ---
+
+# 1. beta.png
+plt.figure()
+plt.plot(f_ext_GHz[:N_m], beta_esprit, label="Forward")
+plt.xlim([8, 12])
+plt.xlabel("Frequency [GHz]")
+plt.ylabel(r"Phase constant $\beta$ [rad/m]")
+plt.title("Dispersion relation")
+plt.savefig("fig1_dispersion.png", dpi=300)
+plt.close()
+
+# 2. alpha.png
+plt.figure()
+plt.plot(f_ext_GHz[:N_m], alpha_esprit, label="Forward")
+plt.xlim([8.5, 12])
+plt.xlabel("Frequency [GHz]")
+plt.ylabel(r"Attenuation rate $\alpha$ [Np/m]")
+plt.title("Attenuation Spectrum")
+plt.savefig("fig2_attenuation.png", dpi=300)
+plt.close()
+
+# 3. ey.png
+plt.figure()
+plt.plot(ZZ_axis * 1000, np.abs(ey_extff_rec[:, Nf_ey]))
+plt.ylabel("|Ey| [V/m]")
+plt.xlabel("z position [mm]")
+plt.title(f"{f_ext_GHz[Nf_ey]:.2f} GHz")
+plt.savefig("fig3_ey.png", dpi=300)
+plt.close()
+
+# --- Interpolation and Export ---
+f_start, f_end = 8.5e9, 12e9
+Nf_start = int(round(f_start / f_ext_step))
+Nf_end = int(round(f_end / f_ext_step))
+
+f_interp = np.linspace(f_ext_GHz[Nf_start], f_ext_GHz[Nf_end], 150)
+f_source = f_ext_GHz[Nf_start : Nf_end + 1]
+alpha_source = alpha_m_db[Nf_start : Nf_end + 1]
+
+spline_interp = interp1d(f_source, alpha_source, kind="cubic", fill_value="extrapolate")
+alpha_total = 4.9 * spline_interp(f_interp)
+
+# 4. alpha_total.png
+plt.figure()
+plt.plot(f_interp, alpha_total, "k")
+plt.xlim([8.5, 12])
+plt.ylim([0, 90])
+plt.xlabel("Frequency [GHz]")
+plt.ylabel("Attenuation over 49mm [dB]")
+plt.title("Attenuation Spectrum")
+plt.savefig("fig4_atten_total.png", dpi=300)
+plt.close()
+
+alpha_tot_data = np.column_stack((f_interp, alpha_total))
+np.savetxt("alpha_total.txt", alpha_tot_data, header="Frequency Attenuation", fmt="%.6f", comments="")

Examples/Artemis_Tutorial_Cases/batch_convert.py

@@ -0,0 +1,43 @@
+import glob
+import os
+import re
+import numpy as np
+import yt
+
+# Directory containing your plotfile folders (e.g., plt00001000, plt00002000, etc.)
+input_dir = "./diags_LLG"
+output_dir = "./raw_data_LLG"
+
+os.makedirs(output_dir, exist_ok=True)
+
+# Find and sort all plotfile folders numerically by timestep
+plt_folders = sorted(
+    [d for d in glob.glob(os.path.join(input_dir, "plt*")) if os.path.isdir(d)],
+    key=lambda x: int(re.search(r'plt(\d+)', x).group(1))
+)
+
+print(f"Found {len(plt_folders)} plotfiles to process.")
+
+for plt_path in plt_folders:
+    # Extract integer timestep from folder name (e.g., 'plt00001000' -> 1000)
+    step_num = int(re.search(r'plt(\d+)', plt_path).group(1))
+
+    # Load dataset and extract covering grid
+    ds = yt.load(plt_path)
+    ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
+
+    # Extract 3D Ey array and slice 1D profile at transverse cell (0, 0)
+    E_array = ad0['Ey'].to_ndarray()
+    ey_1d = E_array[0, 0, :]
+
+    # Generate 0-based k indices
+    k_indices = np.arange(len(ey_1d))
+
+    # Format 2-column data matrix: [k_index, Ey]
+    output_data = np.column_stack((k_indices, ey_1d))
+
+    # Save output text file matching your naming scheme (e.g., 1000.txt)
+    output_filename = os.path.join(output_dir, f"{step_num}.txt")
+    np.savetxt(output_filename, output_data, fmt="%d  %12.8e")
+
+print(f"Batch conversion complete! Output saved to '{output_dir}'.")

Examples/Waveguide/inputs_3d_LLG_filter → Examples/Artemis_Tutorial_Cases/inputs_3d_LLG_filter

@@ -101,10 +101,10 @@ warpx.My_external_grid_function(x,y,z)= "1.3926e5 * (x<=thickness-width/2) + 0 *
 warpx.Mz_external_grid_function(x,y,z) = "0."
 
 #Diagnostics
-diagnostics.diags_names = zline_BernadoFilter0111_
+diagnostics.diags_names = plt
 
-zline_BernadoFilter0111_.intervals = 2
-zline_BernadoFilter0111_.diag_lo = 0.0 0.0 -250.e-3
-zline_BernadoFilter0111_.diag_hi = 0.0 0.0  250.e-3
-zline_BernadoFilter0111_.diag_type = Full
-zline_BernadoFilter0111_.fields_to_plot = Ey Hx Hz
+plt.intervals = 20
+plt.diag_lo = 0.0 0.0 -250.e-3
+plt.diag_hi = 0.0 0.0  250.e-3
+plt.diag_type = Full
+plt.fields_to_plot = Ey Hx Hz

Examples/Waveguide/inputs_3d_empty_X_band → Examples/Artemis_Tutorial_Cases/inputs_3d_empty_X_band

@@ -8,18 +8,18 @@
 ####### GENERAL PARAMETERS ######
 #################################
 max_step = 1000000
-amr.n_cell = 128 4 128 # number of cells spanning the domain in each coordinate direction at level 0
+amr.n_cell = Nx Ny Nz # number of cells spanning the domain in each coordinate direction at level 0
 amr.max_grid_size = 64 # maximum size of each AMReX box, used to decompose the domain
 amr.blocking_factor = 4 # only meaningful for AMR
 geometry.dims = 3
-boundary.field_lo = pec pec pml  # PEC on side walls; PML at -z end
-boundary.field_hi = pec pec pec  # PEC on side walls; PEC at +z end
+boundary.field_lo = pec pec pml # PEC on side walls; PML at -z end
+boundary.field_hi = pec pec pec # PEC on side walls; PEC at +z end
 
 # waveguide width 14.95mm, height is 11.43mm, and length is 500mm
 # 14.95mm to make the wave guide work on fundamental mode at 10.5GHz
 # 11.43mm < 14.95mm so that the fundamental mode is TE10
-geometry.prob_lo = -7.475e-3 -5.715e-3 -250.0e-3 # must be consistent with my_constants.length and .diag_hi/lo
-geometry.prob_hi =  7.475e-3  5.715e-3  250.0e-3
+geometry.prob_lo = -width/2 -height/2 -length/2 # must be consistent with my_constants.length and .diag_hi/lo
+geometry.prob_hi = width/2 height/2 length/2
 
 amr.max_level = 0
 
@@ -28,14 +28,23 @@ my_constants.c = 299792458.
 my_constants.width = 14.95e-3 # waveguide width is 14.95mm
 my_constants.height = 11.43e-3 # waveguide height is 10.16mm
 my_constants.length = 500.0e-3 # waveguide length is 400mm
-my_constants.rjx = 1.0e-4 # the x dimension of current source cross-section
-my_constants.rjz = 10.0e-4 # the z dimension of current source cross-section; should be just larger than 2*dz
 my_constants.wavelength = 0.0286 # frequency is 10.5 GHz
 my_constants.TP = 9.5238e-11 # Gaussian pulse width, 1 x time period of excitation
 my_constants.flag_none = 0 # no source flag
 my_constants.flag_hs = 1 # hard source flag
 my_constants.flag_ss = 2 # soft source flag
 
+my_constants.Nx = 128
+my_constants.Ny = 4
+my_constants.Nz = 128
+
+my_constants.dx = width/Nx
+my_constants.dy = height/Ny
+my_constants.dz = length/Nz
+
+my_constants.rjx = dx/100
+my_constants.rjz = dz/100
+
 #################################
 ############ NUMERICS ###########
 #################################
@@ -58,19 +67,27 @@ macroscopic.mu_function(x,y,z) = "1.25663706212e-06"
 #################################
 
 warpx.B_excitation_on_grid_style = "parse_B_excitation_grid_function"
-warpx.Bx_excitation_grid_function(x,y,z,t) = "1.25663706212e-06 * 2.5e-5 * (exp(-(t-3*TP)**2/(2*TP**2))*cos(2*pi*c/wavelength*t)) * cos(x/(width/2)*(pi/2)) * (z > - rjz/2 + length/2)"
+warpx.Bx_excitation_grid_function(x,y,z,t) = "1.25663706212e-06 * 2.5e-5 * (exp(-(t-3*TP)**2/(2*TP**2))*cos(2*pi*c/wavelength*t)) * cos(x/(width/2)*(pi/2))"
 # current line source Jy = Bx at the +z end of waveguide; spanning over entire y dimension
 warpx.By_excitation_grid_function(x,y,z,t) = "0.0"
 warpx.Bz_excitation_grid_function(x,y,z,t) = "0.0"
-warpx.Bx_excitation_flag_function(x,y,z) = "flag_ss * (z > - rjz/2 + length/2) + flag_none * (z <= - rjz/2 + length/2)"
+warpx.Bx_excitation_flag_function(x,y,z) = "flag_ss * (z > - rjz + length/2 - dz) + flag_none * (z <= - rjz + length/2 - dz)"
 warpx.By_excitation_flag_function(x,y,z) = "flag_none"
 warpx.Bz_excitation_flag_function(x,y,z) = "flag_none"
 
 #Diagnostics
-diagnostics.diags_names = zline_BernadoFilter0126_
+diagnostics.diags_names = plt chk
+
+plt.intervals = 100
+plt.diag_lo = 0.0 0.0 -length/2
+plt.diag_hi = 0.0 0.0 length/2
+plt.diag_type = Full
+plt.fields_to_plot = Ey Bx
+
+# checkpoint
+chk.intervals = 250000
+chk.diag_type = Full
+chk.format = checkpoint
+chk.file_min_digits = 7
 
-zline_BernadoFilter0126_.intervals = 40
-zline_BernadoFilter0126_.diag_lo = 0.0 0.0 -250.e-3
-zline_BernadoFilter0126_.diag_hi = 0.0 0.0  250.e-3
-zline_BernadoFilter0126_.diag_type = Full
-zline_BernadoFilter0126_.fields_to_plot = Ey Bx
+# amr.restart = diags/chk0500000