From 12b7c03a52533dc3eb02f5f2e1bd38407bf6b565 Mon Sep 17 00:00:00 2001 From: Pralay Vaggu Date: Mon, 15 Sep 2025 12:13:09 -0400 Subject: [PATCH] added ALFs example --- init/ALFs_STEVE/calc_Epar.m | 250 +++++++++++++++++++++ init/ALFs_STEVE/calc_Epar.py | 73 ++++++ init/ALFs_STEVE/config.nml | 73 ++++++ init/ALFs_STEVE/fac_said.m | 45 ++++ init/ALFs_STEVE/fac_said.py | 27 +++ init/ALFs_STEVE/inputs/dat2h5_conversion.m | 11 + init/ALFs_STEVE/inputs/setup_grid.json | 16 ++ init/ALFs_STEVE/inputs/simsize.h5 | Bin 0 -> 6156 bytes init/ALFs_STEVE/pot_said.py | 54 +++++ 9 files changed, 549 insertions(+) create mode 100644 init/ALFs_STEVE/calc_Epar.m create mode 100644 init/ALFs_STEVE/calc_Epar.py create mode 100644 init/ALFs_STEVE/config.nml create mode 100644 init/ALFs_STEVE/fac_said.m create mode 100644 init/ALFs_STEVE/fac_said.py create mode 100644 init/ALFs_STEVE/inputs/dat2h5_conversion.m create mode 100644 init/ALFs_STEVE/inputs/setup_grid.json create mode 100644 init/ALFs_STEVE/inputs/simsize.h5 create mode 100644 init/ALFs_STEVE/pot_said.py diff --git a/init/ALFs_STEVE/calc_Epar.m b/init/ALFs_STEVE/calc_Epar.m new file mode 100644 index 0000000..2b6eb05 --- /dev/null +++ b/init/ALFs_STEVE/calc_Epar.m @@ -0,0 +1,250 @@ +% Set the GEMINI_ROOT path +setenv('GEMINI_ROOT', '/Users/vaggu/Work/gemini3d/build/msis/'); + +% Read simulation data +iframe = 21; +direc = '/Users/vaggu/Work/GEMINI_projects/ALFs_STEVE/Non-linear-currents/'; +cfg = gemini3d.read.config(direc); +xg = gemini3d.read.grid(direc); +dat = gemini3d.read.frame(direc,"time",cfg.times(iframe)); % Use frame 16 +y = xg.x3(3:end-2); +z = xg.x1(3:end-2); + +% Recompute conductivities over grid +[sigP, sigH, sig0, SigP, SigH, incap, Incap] = ... + gemini3d.gemscr.postprocess.conductivity_reconstruct(cfg.times(iframe), dat, cfg, xg); + +% Have J1, can compute sigma0, so can get E1 +Jz = dat.J1; % Field-aligned current, positive up +sig0 = permute(sig0,[1 3 2]); +sigH = permute(sigH,[1 3 2]); +sigP = permute(sigP,[1 3 2]); +Ez = Jz ./ sig0; + +% compute E_perp in E-W and N-S +%using conductivity tensor +% Inputs (array-sized): J2 (E–W), J3 (N–S), sigP, sigH +dnm = sigP.^2 + sigH.^2 ; +E2 = (-sigH .* dat.J3 + sigP .* dat.J2) ./ dnm; % E–W electric field +E3 = ( sigP .* dat.J3 + sigH .* dat.J2) ./ dnm; % N–S electric field + +%% plot E, sigma, and J (both perp and par) +figure('Units', 'Inches', 'Position', [0, 0, 16, 16], 'PaperPositionMode', 'auto'); + +%electric fields +subplot(431) +Ezplot = reshape(Ez, [length(z), length(y)]); +pcolor(y./1e3, z./1e3, Ezplot); +shading interp; +% xlabel('northward dist. [km]'); +ylabel('altitude [km]'); +ylim([80, 500]); +xlim([-7, 5]); +colorbar; +title('$E_\parallel$ [V/m]', 'Interpreter', 'latex'); +% grid on; +ax = gca; +% ax.GridLineStyle = ':'; +% ax.GridColor = 'k'; +% ax.GridAlpha = 0.5; +% ax.LineWidth = 1.2; +ax.FontSize = 14; + +subplot(432) +Ezplot = reshape(E2, [length(z), length(y)]); +pcolor(y./1e3, z./1e3, Ezplot); +shading interp; +% xlabel('northward dist. [km]'); +% ylabel('altitude [km]'); +ylim([80, 500]); +xlim([-7, 5]); +colorbar; +title('$E_\perp$ (E-W) [V/m]', 'Interpreter', 'latex'); +% grid on; +ax = gca; +% ax.GridLineStyle = ':'; +% ax.GridColor = 'k'; +% ax.GridAlpha = 0.5; +% ax.LineWidth = 1.2; +ax.FontSize = 14; + +subplot(433) +Ezplot = reshape(E3, [length(z), length(y)]); +pcolor(y./1e3, z./1e3, Ezplot); +shading interp; +% xlabel('northward dist. [km]'); +% ylabel('altitude [km]'); +ylim([80, 500]); +xlim([-7, 5]); +colorbar; +title('$E_\perp$ (N-S) [V/m]', 'Interpreter', 'latex'); +% grid on; +ax = gca; +% ax.GridLineStyle = ':'; +% ax.GridColor = 'k'; +% ax.GridAlpha = 0.5; +% ax.LineWidth = 1.2; +ax.FontSize = 14; + +%currents +subplot(434) +J1plot = reshape(dat.J1, [length(z), length(y)]); +pcolor(y./1e3, z./1e3, J1plot); +shading interp; +% xlabel('northward dist. [km]'); +ylabel('altitude [km]'); +ylim([80, 500]); +xlim([-7, 5]); +colorbar; +title('$J_\parallel$ [A/m$^2$]', 'Interpreter', 'latex'); +ax = gca; +% ax.GridLineStyle = ':'; +% ax.GridColor = 'k'; +% ax.GridAlpha = 0.5; +% ax.LineWidth = 1.2; +ax.FontSize = 14; + +subplot(435) +J2plot = reshape(dat.J2, [length(z), length(y)]); +pcolor(y./1e3, z./1e3, J2plot); +shading interp; +% xlabel('northward dist. [km]'); +% ylabel('altitude [km]'); +ylim([80, 500]); +xlim([-7, 5]); +colorbar; +title('$J_\perp$ (E-W) [A/m$^2$]', 'Interpreter', 'latex'); +ax = gca; +% ax.GridLineStyle = ':'; +% ax.GridColor = 'k'; +% ax.GridAlpha = 0.5; +% ax.LineWidth = 1.2; +ax.FontSize = 14; + +subplot(436) +J3plot = reshape(dat.J3, [length(z), length(y)]); +pcolor(y./1e3, z./1e3, J3plot); +shading interp; +% xlabel('northward dist. [km]'); +% ylabel('altitude [km]'); +ylim([80, 500]); +xlim([-7, 5]); +colorbar; +title('$J_\perp$ (N-S) [A/m$^2$]', 'Interpreter', 'latex'); +ax = gca; +% ax.GridLineStyle = ':'; +% ax.GridColor = 'k'; +% ax.GridAlpha = 0.5; +% ax.LineWidth = 1.2; +ax.FontSize = 14; + +%conductivities +subplot(437) +sig0plot = reshape(sig0, [length(z), length(y)]); +pcolor(y./1e3, z./1e3, sig0plot); +shading interp; +% xlabel('northward dist. [km]'); +ylabel('altitude [km]'); +ylim([80, 500]); +xlim([-7, 5]); +colorbar; +title('$Sig0_\parallel$ [S/m]', 'Interpreter', 'latex'); +% grid on; +ax = gca; +% ax.GridLineStyle = ':'; +% ax.GridColor = 'k'; +% ax.GridAlpha = 0.5; +% ax.LineWidth = 1.2; +ax.FontSize = 14; + +subplot(438) +sigHplot = reshape(sigH, [length(z), length(y)]); +pcolor(y./1e3, z./1e3, sigHplot); +shading interp; +% xlabel('northward dist. [km]'); +% ylabel('altitude [km]'); +ylim([80, 500]); +xlim([-7, 5]); +colorbar; +title('$SigH_\perp$ [S/m]', 'Interpreter', 'latex'); +% grid on; +ax = gca; +% ax.GridLineStyle = ':'; +% ax.GridColor = 'k'; +% ax.GridAlpha = 0.5; +% ax.LineWidth = 1.2; +ax.FontSize = 14; + +subplot(439) +sigPplot = reshape(sigP, [length(z), length(y)]); +pcolor(y./1e3, z./1e3, sigPplot); +shading interp; +% xlabel('northward dist. [km]'); +% ylabel('altitude [km]'); +ylim([80, 500]); +xlim([-7, 5]); +colorbar; +title('$SigP_\perp$ [S/m]', 'Interpreter', 'latex'); +% grid on; +ax = gca; +% ax.GridLineStyle = ':'; +% ax.GridColor = 'k'; +% ax.GridAlpha = 0.5; +% ax.LineWidth = 1.2; +ax.FontSize = 14; + +%ion-electron mobilities +subplot(4,3,10) +v1plot = reshape(dat.v1, [length(z), length(y)]); +pcolor(y./1e3, z./1e3, v1plot); +shading interp; +xlabel('northward dist. [km]'); +ylabel('altitude [km]'); +ylim([80, 500]); +xlim([-7, 5]); +colorbar; +title('$v_\parallel$ [m/s]', 'Interpreter', 'latex'); +% grid on; +ax = gca; +% ax.GridLineStyle = ':'; +% ax.GridColor = 'k'; +% ax.GridAlpha = 0.5; +% ax.LineWidth = 1.2; +ax.FontSize = 14; + +subplot(4,3,11) +pcolor(y./1e3, z./1e3, dat.v2); +shading interp; +xlabel('northward dist. [km]'); +% ylabel('altitude [km]'); +ylim([80, 500]); +xlim([-7, 5]); +colorbar; +title('$v_\perp$ (E-W) [m/s]', 'Interpreter', 'latex'); +% grid on; +ax = gca; +% ax.GridLineStyle = ':'; +% ax.GridColor = 'k'; +% ax.GridAlpha = 0.5; +% ax.LineWidth = 1.2; +ax.FontSize = 14; + +subplot(4,3,12) +pcolor(y./1e3, z./1e3, dat.v3); +shading interp; +xlabel('northward dist. [km]'); +% ylabel('altitude [km]'); +ylim([80, 500]); +xlim([-7, 5]); +colorbar; +title('$v_\perp$ (N-S) [m/s]', 'Interpreter', 'latex'); +% grid on; +ax = gca; +% ax.GridLineStyle = ':'; +% ax.GridColor = 'k'; +% ax.GridAlpha = 0.5; +% ax.LineWidth = 1.2; +ax.FontSize = 14; +p0 = 'E_J_V_sigma.png'; +fullpath = fullfile(direc,p0); +saveas(gcf, fullpath); \ No newline at end of file diff --git a/init/ALFs_STEVE/calc_Epar.py b/init/ALFs_STEVE/calc_Epar.py new file mode 100644 index 0000000..61fd2e6 --- /dev/null +++ b/init/ALFs_STEVE/calc_Epar.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Jul 17 14:32:28 2024 + +Because GEMINI2D does not output Eparallel or full potential (i.e. varying along B) + by default we need to do a few calculations to recover these from the data. + +NOTES: + - add potential drop calculations + - add mean free path calculation, total potential drop is irrelevant if all the + energy is lost soon after acceleration + - temperature may play a role in mfp, i.e. both in thermal velocity and also in + temp.-dependence of collisions (need high energy corrections?) + +@author: zettergm +""" + +import gemini3d.read +import gemini3d.conductivity +import matplotlib.pyplot as plt +import os +import numpy as np + +############################################################################## +# tell pygemini where to find msis +os.environ["GEMINI_ROOT"]="/Users/zettergm/Projects/gemini3d/build/msis/" + +# read simulations data +iframe=16 +direc="~/simulations/sdcard/STEVE2D_dist_test/" +cfg=gemini3d.read.config(direc) +xg=gemini3d.read.grid(direc) +dat=gemini3d.read.frame(direc,cfg["time"][iframe]) # use frame 50 +y=xg["x3"][2:-2] +z=xg["x1"][2:-2] + +# recompute conductivities over grid +sigP, sigH, sig0, SigP, SigH, incap, Incap = ( + gemini3d.conductivity.conductivity_reconstruct(cfg["time"][iframe], dat, cfg, xg) ) + +# have J1, can compute sigma0, so can get E1 +Jz=dat["J1"] # field-aligned current, positive up +Ez=Jz/sig0 +############################################################################## + +plt.figure(dpi=200) +Ezplot=np.reshape(np.array(Ez), [z.size,y.size]) +plt.pcolormesh(y,z,Ezplot,shading="gouraud") +plt.xlabel("y") +plt.ylabel("z") +plt.ylim((50e3, 400e3)) +plt.colorbar() +plt.title("$E_\parallel$ (V/m)") + +plt.figure(dpi=200) +neplot=np.reshape(np.array(dat["ne"]), [z.size,y.size]) +plt.pcolormesh(y,z,np.log10(neplot),shading="gouraud") +plt.xlabel("y") +plt.ylabel("z") +plt.ylim((50e3, 400e3)) +plt.clim(4,11.5) +plt.colorbar() +plt.title("$n_e$ (V/m)") + +plt.figure(dpi=200) +J1plot=np.reshape(np.array(dat["J1"]), [z.size,y.size]) +plt.pcolormesh(y,z,J1plot,shading="gouraud") +plt.xlabel("y") +plt.ylabel("z") +plt.ylim((50e3, 400e3)) +plt.colorbar() +plt.title("$J_\parallel$ (V/m)") \ No newline at end of file diff --git a/init/ALFs_STEVE/config.nml b/init/ALFs_STEVE/config.nml new file mode 100644 index 0000000..d154a99 --- /dev/null +++ b/init/ALFs_STEVE/config.nml @@ -0,0 +1,73 @@ +&base +ymd = 2008,3,26 ! year, month, day +UTsec0 = 12600.0 ! start time in UT seconds +tdur = 300.0 ! duration of simulation in seconds +dtout = 15.0 ! how often to do output +activ = 150.0,150.0,50.0 ! f107a,f107,Ap +tcfl = 0.9 ! target cfl number +Teinf = 1500.0 ! exospheric electron temperature +/ + +&setup +glat = 65 +glon = 213 +xdist = 1500e3 ! eastward distance (meters) +ydist = 20e3 ! northward distance (meters) +lxp = 1 +lyp = 1032 +Bincl = 90 +nmf = 5e11 +nme = 2e11 +!alt_min = 80e3 +alt_min = 80e3 +alt_max = 950e3 +!alt_scale = 10e3, 9.9e3, 500e3, 150e3 ! altitude grid scales (meters) +alt_scale = 10e3, 9.9e3, 500e3, 150e3 ! altitude grid scales (meters) +eq_dir = "~/Work/LWS_ARCS/GEMINI/arcs_eq" +!Efield_latwidth = 0.0125 +Efield_latwidth = 0.0250 +!Jtarg = 100e-6 +Jtarg = 25e-6 +Jtarg_function = "fac_said" +Efield_llat = 1024 +/ + +&flags +potsolve = 1 ! solve electrodynamics: 0 - no; 1 - electrostatic; 2 - inductive +flagoutput = 2 +/ + +&files +indat_size = 'inputs/simsize.h5' +indat_grid = 'inputs/simgrid.h5' +indat_file = 'inputs/initial_conditions.h5' +/ + +&efield +dtE0 = 60.0 +E0_dir = 'inputs/Efield/' +/ + +!&precip_BG +!PhiWBG=5e-1 ! total energy flux (mW/m^2) +!W0BG=1e3 ! characteristic energy (eV) +!/ + +&precip_BG +PhiWBG=100 ! total energy flux (mW/m^2) +!PhiWBG=0.5 ! total energy flux (mW/m^2) +W0BG=3000.0 ! characteristic energy (eV) +/ + +&flags +flagFBI = 2 ! Farley-Buneman Instability: 0 - no; 1 - Only anomalous heating; 2 - Anomalous heating + nonlinear current +/ + +! (optional - default value of 2) +&diffusion +diffsolvetype = 1 ! type of diffusion solver to use: 1 - backward Euler; 2 - TRBDF2 +/ + +&milestone +mcadence = 4 +/ diff --git a/init/ALFs_STEVE/fac_said.m b/init/ALFs_STEVE/fac_said.m new file mode 100644 index 0000000..27ccd5c --- /dev/null +++ b/init/ALFs_STEVE/fac_said.m @@ -0,0 +1,45 @@ +function E = fac_said(E, Nt, gridflag, flagdip) + % Function for 3D simulations, calculates FAC up/down 0.5 degree FWHM + % Inputs: + % E - Structure or Dataset (similar to xarray.Dataset) + % gridflag - integer (used to determine direction) + % flagdip - boolean flag + + % Check if E.mlon or E.mlat has only one value (similar to xarray dimension checks) + % if length(E.mlon) == 1 || length(E.mlat) == 1 + % error('This function is for 3D simulations only.'); + % end + + % Uniform in longitude + shapelon = 1; + + config = "one-sided"; + % Calculate shape in latitude (similar to NumPy operations) + if config == "one-sided" + shapelat = exp(-( (E.mlat - E.mlatmean - 1.5 * E.mlatsig).^2 ) / (2 * E.mlatsig^2)) ... + - exp(-( (E.mlat - E.mlatmean + 1.5 * E.mlatsig).^2 ) / (2 * E.mlatsig^2)); + elseif config == "two-sided" + shapelat = 1/2*exp(-( (E.mlat - E.mlatmean - 1.5 * E.mlatsig).^2 ) / (2 * E.mlatsig^2)) ... + - exp(-( (E.mlat - E.mlatmean + 1.5 * E.mlatsig).^2 ) / (2 * E.mlatsig^2))... + + 1/2*exp(-( (E.mlat - E.mlatmean + 4.5 * E.mlatsig).^2 ) / (2 * E.mlatsig^2)); + end + + + % Loop over times in E.time, starting from the 3rd element (MATLAB indices start at 1) + for t_idx = 3:length(E.times) + t = E.times(t_idx); + + % Set flagdiriich to zero + E.flagdirich(t_idx) = 2; + + % Set key based on gridflag value + if gridflag == 1 + k = 'Vminx1it'; + else + k = 'Vmaxx1it'; + end + + % Update E.(k) based on the calculations + E.(k)(:,:,t_idx) = E.Jtarg * shapelon * shapelat; + end +end \ No newline at end of file diff --git a/init/ALFs_STEVE/fac_said.py b/init/ALFs_STEVE/fac_said.py new file mode 100644 index 0000000..276bbfb --- /dev/null +++ b/init/ALFs_STEVE/fac_said.py @@ -0,0 +1,27 @@ +#STEVE_py_func +import xarray +import numpy as np + + +def fac_said(E: xarray.Dataset, gridflag: int, flagdip: bool) -> xarray.Dataset: + """ + for 3D sim, FAC up/down 0.5 degree FWHM + """ + + # if E.mlon.size == 1 or E.mlat.size == 1: + # raise ValueError("for 3D sims only") + + # uniform in longitude + shapelon = 1 + + shapelat = np.exp( + -((E.mlat - E.mlatmean - 1.5 * E.mlatsig) ** 2) / 2 / E.mlatsig**2 + ) - np.exp(-((E.mlat - E.mlatmean + 1.5 * E.mlatsig) ** 2) / 2 / E.mlatsig**2) + + for t in E.time[2:]: + E["flagdirich"].loc[t] = 0 + + k = "Vminx1it" if gridflag == 1 else "Vmaxx1it" + E[k].loc[t] = E.Jtarg * shapelon * shapelat + + return E diff --git a/init/ALFs_STEVE/inputs/dat2h5_conversion.m b/init/ALFs_STEVE/inputs/dat2h5_conversion.m new file mode 100644 index 0000000..bd7dd53 --- /dev/null +++ b/init/ALFs_STEVE/inputs/dat2h5_conversion.m @@ -0,0 +1,11 @@ +close all;clear all;clc; + +filename_dat = "initial_conditions.dat"; +fileID = fopen(filename_dat,'r'); + +data = fread(fileID,'double'); +fclose(fileID); + +filename_h5 = "initial_conditions.h5"; +h5create(filename_h5,'/dataset1',size(data)); +h5write(filename_h5,'/dataset1',data); diff --git a/init/ALFs_STEVE/inputs/setup_grid.json b/init/ALFs_STEVE/inputs/setup_grid.json new file mode 100644 index 0000000..fd83276 --- /dev/null +++ b/init/ALFs_STEVE/inputs/setup_grid.json @@ -0,0 +1,16 @@ +{ + "matlab": { + "arch": "maci64", + "version": "9.14.0.2674353 (R2023a) Update 7" + }, + "git": { + "version": "git version 2.39.3 (Apple Git-146)", + "remote": "https://github.com/gemini3d/mat_gemini", + "branch": "main", + "commit": "v5.6.0-101-g089c338", + "porcelain": "false" + }, + "eq": { + "eq_dir": "/Users/vaggu/Work/Ph.D./LWS_ARCS/GEMINI/arcs_eq" + } +} \ No newline at end of file diff --git a/init/ALFs_STEVE/inputs/simsize.h5 b/init/ALFs_STEVE/inputs/simsize.h5 new file mode 100644 index 0000000000000000000000000000000000000000..9ca9e9bec9b5966f598a2bc4c11fbecc659d3fbb GIT binary patch literal 6156 zcmeHKIS#@w5L_n#SrJG;g{blbj)oSAASysWK}kuG2jmIR@c6JAfj1P{hX~!N3I12V z^)gM#crrgD=%^cci+Q+S@IMrorlW}EIb_E#QBl>! zug`3R>L xarray.Dataset: + """ + synthesize a feature + """ + + if E.Etarg > 1: + logging.warning(f"Etarg units V/m -- is {E['Etarg']} V/m realistic?") + + # NOTE: h2, h3 have ghost cells, so we use lx1 instead of -1 to index + # pk is a scalar. + # north-south + S = E.Etarg * E.sigx3 * xg["h3"][lx1, 0, lx3 // 2] * np.sqrt(np.pi) / 2 + mlatmean = np.mean(E.mlat) + + taper = erf((E.mlat - mlatmean - E.mlatoffset) / E.mlatsig).data[None, :] + taper = taper - erf((E.mlat - mlatmean) / E.mlatoffset).data[None, :] + 1 + taper = taper + erf((E.mlat - mlatmean + E.mlatoffset) / E.mlatsig).data[None, :] + taper = -1 * taper + + assert S.ndim == 0, "S is a scalar" + + for t in E.time: + E["flagdirich"].loc[t] = 1 + + if gridflag == 1: + E["Vminx1it"].loc[t] = S * taper + else: + E["Vmaxx1it"].loc[t] = S * taper + + return E