Skip to content

Commit

Permalink
radiation: add some basic dipole / quadrupole radiation
Browse files Browse the repository at this point in the history
This really could use some more work in cleaning it up, etc, but
it does pretty much do what it's supposed to do.
  • Loading branch information
germasch committed Nov 7, 2024
1 parent afe118f commit 0dc8cfb
Show file tree
Hide file tree
Showing 2 changed files with 252 additions and 10 deletions.
194 changes: 194 additions & 0 deletions python/psc_radiation.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,194 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib as mpl\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.animation import FuncAnimation\n",
"import xarray as xr\n",
"import numpy as np\n",
"#import pscpy\n",
"\n",
"%matplotlib ipympl \n",
"#plt.rcParams['figure.figsize'] = [16, 10]"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# dipole far field\n",
"dir = \"/Users/kai/src/psc/build-arm64\"\n",
"steps = range(0, 400, 5)\n",
"vmax = .0001\n",
"\n",
"# orig dipole near field\n",
"dir = \"/Users/kai/src/psc/build-arm64-2\"\n",
"steps = range(0, 800, 5)\n",
"vmax = .01\n",
"\n",
"# quadrupole\n",
"dir = \"/Users/kai/src/psc/build-arm64-2\"\n",
"steps = range(0, 1000, 20)\n",
"vmax = .00005\n",
"\n",
"def open_step(step):\n",
" return xr.open_dataset(f\"{dir}/pfd.{step:09d}.bp\", #engine='pscadios2',\n",
" species_names=['e', 'i'])\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"datas = {}\n",
"for step in steps:\n",
" ds = open_step(step)\n",
" datas[step] = ds.ez_ec.sel(y=0.).T, float(ds.time)\n",
"#datas"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def get_data(step):\n",
" return datas[step]\n",
"\n",
"steps = list(steps)\n",
"\n",
"fig, ax = plt.subplots()\n",
"step = steps[0]\n",
"fld, time = get_data(step)\n",
"cax = fld.plot(vmin=-vmax, vmax=vmax, cmap='coolwarm')\n",
"ax.set_title(f\"step {step} time {time:6.2f}\")\n",
"ax.set_aspect(1.)\n",
"\n",
"def animate(step):\n",
" fld, time = get_data(step)\n",
" cax.set_array(fld.values.flatten())\n",
" ax.set_title(f\"step {step} time {time:6.2f}\")\n",
"\n",
"ani = FuncAnimation(\n",
" fig, # figure\n",
" animate, # name of the function above\n",
" frames=steps[1:], # Could also be iterable or list\n",
" interval=200, # ms between frames\n",
" blit=False\n",
")\n",
"\n",
"plt.show()\n",
"#ani"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def plot_fields(fldnames, fld_kwargs=None):\n",
" fig, axs = plt.subplots(1, len(fldnames))\n",
" if len(fldnames) == 1: axs = [axs]\n",
" for i, fldname in enumerate(fldnames):\n",
" fld = ds[fldname].sel(y=0)\n",
" if fld_kwargs:\n",
" kwargs = fld_kwargs[i]\n",
" else:\n",
" kwargs = {}\n",
" fld.plot(ax=axs[i], **kwargs)\n",
" axs[i].set_aspect('equal')\n",
" plt.tight_layout()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"plot_fields(['ex_ec', 'ey_ec', 'ez_ec'])\n",
"# fld_kwargs=[{\"vmin\": -.0065}, {\"vmin\": -.02}, {}])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def format_radians_label(float_in):\n",
" # Converts a float value in radians into a\n",
" # string representation of that float\n",
" string_out = str(float_in / (np.pi))+\"π\"\n",
" \n",
" return string_out\n",
"\n",
"def convert_polar_xticks_to_radians(ax):\n",
" # Converts x-tick labels from degrees to radians\n",
" \n",
" # Get the x-tick positions (returns in radians)\n",
" label_positions = ax.get_xticks()\n",
" \n",
" # Convert to a list since we want to change the type of the elements\n",
" labels = list(label_positions)\n",
" \n",
" # Format each label (edit this function however you'd like)\n",
" labels = [format_radians_label(label) for label in labels]\n",
" \n",
" ax.set_xticklabels(labels)\n",
" \n",
"theta = np.linspace(-np.pi, np.pi, 100)\n",
"\n",
"fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})\n",
"ax.plot(theta, np.sin(theta)**2)\n",
"ax.set_rticks([0.25, 0.5, 0.75, 1])\n",
"ax.set_theta_zero_location(\"N\")\n",
"#convert_polar_xticks_to_radians(ax)\n",
"\n",
"ax.set_title(\"Radiated dipole power\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": ".venv",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.5"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
68 changes: 58 additions & 10 deletions src/psc_radiation.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,9 @@ struct PscFlatfoilParams
double k = 2 * 2. * M_PI / 10.;
double amplitude_s = 1.;
double amplitude_p = 0.;

double omega = 1.;
double d = .1;
};

// ======================================================================
Expand Down Expand Up @@ -84,9 +87,9 @@ using Moment_n = typename Moment_n_Selector<Mparticles, Dim>::type;
void setupParameters()
{
// -- set some generic PSC parameters
psc_params.nmax = 10000001; // 5001;
psc_params.nmax = 1001; // 5001;
psc_params.cfl = 0.75;
psc_params.write_checkpoint_every_step = 1000;
// psc_params.write_checkpoint_every_step = 1000;
psc_params.stats_every = 1;

// -- start from checkpoint:
Expand All @@ -102,7 +105,7 @@ void setupParameters()
// -- Set some parameters specific to this case
g.theta = 0;
g.k = 2 * 2. * M_PI / 10.;
g.amplitude_s = 1.;
g.amplitude_s = 0.;
g.amplitude_p = 0.;
}

Expand All @@ -117,9 +120,16 @@ void setupParameters()
Grid_t* setupGrid()
{
// --- setup domain
Grid_t::Real3 LL = {10, 10., 10.}; // domain size (in d_e)
Int3 gdims = {50, 50, 50}; // global number of grid points
Int3 np = {5, 5, 5}; // division into patches
// far field
// Grid_t::Real3 LL = {80., 80., 80.}; // domain size (in d_e)
// Int3 gdims = {400, 400, 400}; // global number of grid points
// near field
// Grid_t::Real3 LL = {5., 5., 5.}; // domain size (in d_e)
// Int3 gdims = {200, 200, 200}; // global number of grid points
// quadrupole far field
Grid_t::Real3 LL = {100., 100., 100.}; // domain size (in d_e)
Int3 gdims = {400, 400, 400}; // global number of grid points
Int3 np = {2, 2, 2}; // division into patches

Grid_t::Domain domain{gdims, LL, -.5 * LL, np};

Expand Down Expand Up @@ -232,6 +242,43 @@ void run()
psc_params.marder_interval = -1;
Marder marder(grid, marder_diffusion, marder_loop, marder_dump);

auto lf_ext_current = [&](const Grid_t& grid, MfieldsState& mflds) {
double time = grid.timestep() * grid.dt;
auto& gdims = grid.domain.gdims;
for (int p = 0; p < mflds.n_patches(); ++p) {
auto& patch = grid.patches[p];
auto F = make_Fields3d<Dim>(mflds[p]);
grid.Foreach_3d(0, 0, [&](int i, int j, int k) {
Int3 index{i, j, k};
auto crd_ec_z = Centering::getPos(patch, index, Centering::EC, 2);
#if 0
double r =
std::sqrt(sqr(crd_ec_z[0]) + sqr(crd_ec_z[1]) + sqr(crd_ec_z[2]));
if (r < g.d) {
F(JZI, i, j, k) += cos(g.omega * time);
}
#else
Int3 ii = Int3{i, j, k} + patch.off;
if (0) { // dipole
if (ii[0] == gdims[0] / 2 && ii[1] == gdims[1] / 2 &&
ii[2] == gdims[2] / 2 - 1) {
F(JZI, i, j, k) += cos(g.omega * time);
}
} else { // quadrupole
if (ii[0] == gdims[0] / 2 && ii[1] == gdims[1] / 2 &&
ii[2] == gdims[2] / 2 - 1) {
F(JZI, i, j, k) += cos(g.omega * time);
}
if (ii[0] == gdims[0] / 2 && ii[1] == gdims[1] / 2 &&
ii[2] == gdims[2] / 2) {
F(JZI, i, j, k) -= cos(g.omega * time);
}
}
#endif
});
}
};

// ----------------------------------------------------------------------
// Set up output
//
Expand All @@ -240,7 +287,7 @@ void run()
// -- output fields
OutputFieldsItemParams outf_item_params{};
OutputFieldsParams outf_params{};
outf_item_params.pfield.out_interval = 1;
outf_item_params.pfield.out_interval = 20;

outf_params.fields = outf_item_params;
OutputFields<MfieldsState, Mparticles, Dim, Writer> outf{grid, outf_params};
Expand All @@ -259,14 +306,15 @@ void run()

if (read_checkpoint_filename.empty()) {
initializeFields(mflds);
lf_ext_current(*grid_ptr, mflds);
}

// ----------------------------------------------------------------------
// hand off to PscIntegrator to run the simulation

auto psc =
makePscIntegrator<PscConfig>(psc_params, *grid_ptr, mflds, mprts, balance,
collision, checks, marder, diagnostics);
auto psc = makePscIntegrator<PscConfig>(
psc_params, *grid_ptr, mflds, mprts, balance, collision, checks, marder,
diagnostics, injectParticlesNone, lf_ext_current);

MEM_STATS();
psc.integrate();
Expand Down

0 comments on commit 0dc8cfb

Please sign in to comment.