diff --git a/python/psc_radiation.ipynb b/python/psc_radiation.ipynb new file mode 100644 index 000000000..446e6a632 --- /dev/null +++ b/python/psc_radiation.ipynb @@ -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 +} diff --git a/src/psc_radiation.cxx b/src/psc_radiation.cxx index 93eb49394..14bf8a2c9 100644 --- a/src/psc_radiation.cxx +++ b/src/psc_radiation.cxx @@ -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; }; // ====================================================================== @@ -84,9 +87,9 @@ using Moment_n = typename Moment_n_Selector::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: @@ -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.; } @@ -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}; @@ -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(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 // @@ -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 outf{grid, outf_params}; @@ -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(psc_params, *grid_ptr, mflds, mprts, balance, - collision, checks, marder, diagnostics); + auto psc = makePscIntegrator( + psc_params, *grid_ptr, mflds, mprts, balance, collision, checks, marder, + diagnostics, injectParticlesNone, lf_ext_current); MEM_STATS(); psc.integrate();