diff --git a/CMakeLists.txt b/CMakeLists.txt index c22ef9de727..dd0a8729c4e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -241,6 +241,9 @@ endif() target_compile_definitions(WarpX PRIVATE WARPX_PARSER_DEPTH=${WarpX_PARSER_DEPTH}) +target_compile_definitions(WarpX PRIVATE PULSAR) + + # Warnings #################################################################### # diff --git a/Examples/Physics_applications/pulsar/Pulsar_scale_RstarPhysical.ipynb b/Examples/Physics_applications/pulsar/Pulsar_scale_RstarPhysical.ipynb new file mode 100644 index 00000000000..f87d443abb6 --- /dev/null +++ b/Examples/Physics_applications/pulsar/Pulsar_scale_RstarPhysical.ipynb @@ -0,0 +1,1392 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Overview\n", + "\n", + "This a notebook that inspects the results of a WarpX simulation.\n", + "\n", + "# Instruction\n", + "\n", + "Enter the path of the data you wish to visualize below. Then execute the cells one by one, by selecting them with your mouse and typing `Shift + Enter`" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "# Import statements\n", + "import yt ; yt.funcs.mylog.setLevel(50)\n", + "import numpy as np\n", + "import scipy.constants as scc\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib notebook" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#Define the Physical Constants, normalizations, and the scale-down parameter, r_scale, for the pulsar. " + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "######################\n", + "# physical constants #\n", + "######################\n", + "c = 299792458.0 # speed of light \n", + "q_e = 1.602176634e-19 # elementary charge\n", + "me=9.10938356*np.power(10.,-31) # electron mass\n", + "epsilon=8.8541878128*np.power(10.,-12) # permittivity of free space\n", + "mu_o = 4*3.14*1e-7 # permeability\n", + "pi = 3.14159265358979323846\n", + "SolarMass = 2e30\n", + "G = 6.674e-11 # gravitational constant\n", + "\n", + "#############################################################################\n", + "# Parameters for a real Pulsar #\n", + "# (Table 7 of J.Petri's article on Theory of Pulsar Magnetosphere and Wind) #\n", + "#############################################################################\n", + "Omega_real = 6283\n", + "B_real = 7.4E4\n", + "R_real = 12000\n", + "n_real = 6.9e16\n", + "omega_pe_real = (n_real*q_e*q_e/(me*epsilon))**0.5 #plasma frequency\n", + "SkinDepth_real = c/omega_pe_real \n", + "Mstar = 1.4*SolarMass\n", + "Rstar_skinDepth_real = 6e5\n", + "\n", + "##################\n", + "# Normalizations #\n", + "##################\n", + "Rstar_skinDepth = 6e0 # Ratio of radius of star to the skin Depth \n", + " # For a real star, this ratio is 6e5\n", + "exponent = np.arange(0,6,1) \n", + "Factor = np.array(10**exponent)\n", + "Rstar_skinDepth = np.array(6*Factor) \n", + "\n", + "RLC_Rstar = 4 # Ratio of light cylinder (where the particle speed ~ c) to Rstar\n", + " # For a pulsar with 1ms period, this ratio is 4. \n", + " # i.e., This ratio sets the omega value, since Omega*RLC = c\n", + "\n", + "# Choose skindepth as the free parameter for a choice of normalizations given above\n", + "# The skin depth below is computed from the number density for a real pulsar\n", + "# Keeping the SkinDepth constant across all the scales in our scaling study\n", + "#SkinDepth = 0.02\n", + "\n", + "# Since skin depth is maintained across the scales, and RLC/Rstar is also maintained \n", + "# lets define the decrease in the scale by comparing the value of \n", + "# (Rstar/skinDepth)/(Rstar_real/skinDepth_real)\n", + "#r_scale = Rstar_skinDepth/Rstar_skinDepth_real\n", + "\n", + "Rstar = np.ones(6)*12000\n", + "SkinDepth = Rstar/Rstar_skinDepth\n", + "r_scale = Rstar_skinDepth_real/Rstar_skinDepth" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "##################################################################\n", + "# Derive other physical parameters from the above normalizations #\n", + "##################################################################\n", + "\n", + "# 1. Lorentz Boost (dimensionless) #\n", + "# Note that in Alex Chen's paper, gamma_o ~ (1/4)*normalized_values.\n", + "# Instead here we have pi/2, since that leads to the closes gamma_o \n", + "# value for a real 1ms pulsar. \n", + "gamma_o = (pi/2)*(Rstar_skinDepth)**2/RLC_Rstar \n", + "\n", + "\n", + "gamma_real = (pi/2)*6e5**2/RLC_Rstar;\n", + "gamma_scaling = gamma_o/(gamma_real) # This is to see how the gamma value \n", + " # decreases due to decrease in the ratio of R_star to skin depth\n", + "\n", + "\n", + "\n", + "# 3. Light cylinder location (m)\n", + "RLC = Rstar*RLC_Rstar\n", + "# 4. Angular Frequency (rad/s)\n", + "# Omega remains constant\n", + "Omega = c/RLC\n", + "# 5. Period (s)\n", + "# Period remains constant\n", + "Period = 2*3.14/Omega\n", + "# Moment of inertia for a sphere = (2/5)MR^2 (kg.m^2)\n", + "# Note that when the Rstar is scaled by r_scale, \n", + "# Mstar decreases as r_scale^3. Thus Mstar*r_scale**3 is the \n", + "# mass of the scaled down star.\n", + "# I remains constant across all scales\n", + "I = (2/5)*(Mstar)*Rstar**2\n", + "# 6. Rotation induced potential drop from pote to equator (V)\n", + "# Reference: Alex Chen's 2014 paper\n", + "# Scales as 1/r_scale**2\n", + "phi_o = gamma_o * (me/q_e) * c**2\n", + "\n", + "# Braking is the rate of slow-down of the spin. \n", + "# It is not relevant for the scaling. \n", + "# However, 1e-15 is the P_dot for a pulsar with P=1s\n", + "# and the rate of slowdown decreases proportional to the period. \n", + "# So we were to compare two stars with same radius, but different Omega\n", + "# then we can use Braking to determine the magnetic field strength at the surface as follows\n", + "# Bo = (3*mu_o*c*c*c/(32*pi*pi*pi))**0.5 * (I*Braking*Period)**0.5/(Rstar**3) \n", + "# (Reference : Table 7 of Jetri's article)\n", + "# Remains constant across all scales\n", + "Braking = 1e-15*Period\n", + "\n", + "# 7. Magnetic field strength (Tesla)\n", + "# Bo decreases as ~ 1/r_scale**2\n", + "Bo = phi_o/(Rstar*Rstar*Omega)\n", + "\n", + "# 8. Volume charge density (C/m^3) for uniform magnetization insize the star \n", + "# Refer to Table 7 of Petri's article \n", + "# Since Omega increases as r_scale and Bo decreases as r_scale\n", + "# The product of Omega*Bo remains constant across all scales. \n", + "# Thus, rho_e decreases as (1/r_scale**2)\n", + "rho_e = 2*epsilon*Omega*Bo #(Not adding the negative sign here)\n", + "\n", + "# 9. Electron number density (#/m^3)\n", + "# ne decreases as (1/r_scale**2)\n", + "ne = rho_e/q_e\n", + "# 9a. plasma frequency \n", + "# decreases as (1/r_scale)\n", + "omega_pe = (ne*q_e*q_e/(me*epsilon))**0.5 #plasma frequency\n", + "# 10. Magnetic moment (Am^2)\n", + "# decreases as (1/r_scale**2)\n", + "magnetic_moment = Bo*Rstar*Rstar*Rstar*4*pi/(mu_o)\n", + "# 11. E-field (V/m)\n", + "# Efield decreases as (1/r_scale**2)\n", + "E = Omega*Bo*Rstar\n", + "\n", + "\n", + "#########################################\n", + "# How do the energies and forces scale? #\n", + "#########################################\n", + "# 12. Magnetic Energy (J)\n", + "# Magnetic energy decreases as (1/r_scale**4)\n", + "magnetic_energy = (4*pi/3)*Bo**2*Rstar**3/(2*mu_o)\n", + "\n", + "# 13. Gravitational Potential Energy (J)\n", + "# G.pot energy remains constant \n", + "GP_energy = (3/5)*G*Mstar*Mstar/(Rstar)\n", + "\n", + "# 14. Gravitational to Electric force (dimensionless)\n", + "# This ratio increases as r_scale**2\n", + "G_EForce = G*Mstar*me/(Rstar*Rstar*q_e*E)\n", + "\n", + "# 15. Rotational kinetic energy \n", + "# Rot. KE remains constant\n", + "rotational_KE = (1/2)*I*Omega**2\n", + "\n", + "# 16. From (12) and (13), we know the B.energy scales as (1/r_scale**4) and GP energy is constant \n", + "# Thus the ratio of GP_energy and B_energy increases as r_scale**4\n", + "GB_energy_ratio = GP_energy/magnetic_energy\n", + "\n", + "# 17. Rate of change of Omega, or angular acceleration decreases as 1/r_scale**4\n", + "Omega_dot = Bo*Bo*Rstar**6*Omega**3/(I) * (32*pi/(3*mu_o*c**3))* (1/(4*pi*pi))\n", + "\n", + "# 18. Torque decreases as 1/r_scale**4\n", + "Torque = I * Omega_dot\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Input parameters for WarpX for geometry scaled down by r_scale = [1.e-05 1.e-04 1.e-03 1.e-02 1.e-01 1.e+00]\n", + "Lorentz factor, (gamma_o) = [1.41371669e+01 1.41371669e+03 1.41371669e+05 1.41371669e+07\n", + " 1.41371669e+09 1.41371669e+11]\n", + "Rstar (m) [12000. 12000. 12000. 12000. 12000. 12000.]\n", + "Omega (rad/s) [6245.67620833 6245.67620833 6245.67620833 6245.67620833 6245.67620833\n", + " 6245.67620833]\n", + "Bo (Tesla) [8.03230942e-06 8.03230942e-04 8.03230942e-02 8.03230942e+00\n", + " 8.03230942e+02 8.03230942e+04]\n", + "ne (/m^3) [5.54482989e+06 5.54482989e+08 5.54482989e+10 5.54482989e+12\n", + " 5.54482989e+14 5.54482989e+16]\n", + "Size of the domain, (m) [360000. 360000. 360000. 360000. 360000. 360000.]\n", + "Minimum cell size (m) [2.e+03 2.e+02 2.e+01 2.e+00 2.e-01 2.e-02]\n", + "timestep (s) [3.76386776e-06 3.76386776e-07 3.76386776e-08 3.76386776e-09\n", + " 3.76386776e-10 3.76386776e-11]\n", + "Numver of cells assuming unif. grid, [1.8e+02 1.8e+03 1.8e+04 1.8e+05 1.8e+06 1.8e+07]\n" + ] + } + ], + "source": [ + "############################################################\n", + "# Print all the values that will be used as input in WarpX #\n", + "############################################################\n", + "print(\"Input parameters for WarpX for geometry scaled down by r_scale = \", Rstar_skinDepth/Rstar_skinDepth_real)\n", + "print(\"Lorentz factor, (gamma_o) = \", gamma_o)\n", + "print(\"Rstar (m)\", Rstar)\n", + "print(\"Omega (rad/s)\", Omega)\n", + "print(\"Bo (Tesla)\", Bo)\n", + "print(\"ne (/m^3)\", ne)\n", + "print(\"Size of the domain, (m)\", 30*Rstar)\n", + "print(\"Minimum cell size (m)\", SkinDepth)\n", + "print(\"timestep (s)\", 0.5*omega_pe**-1) # Or may be cfl criterion\n", + "print(\"Numver of cells assuming unif. grid, \", 30*Rstar/SkinDepth)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Gravitational force to E-force : [1.22562947e-02 1.22562947e-04 1.22562947e-06 1.22562947e-08\n", + " 1.22562947e-10 1.22562947e-12]\n", + "Gravitational energy to B-energy : [1.40727411e+38 1.40727411e+34 1.40727411e+30 1.40727411e+26\n", + " 1.40727411e+22 1.40727411e+18]\n", + "B energy, (J) [1.85906071e+08 1.85906071e+12 1.85906071e+16 1.85906071e+20\n", + " 1.85906071e+24 1.85906071e+28]\n", + "Gravitational potential energy, (J) [2.616208e+46 2.616208e+46 2.616208e+46 2.616208e+46 2.616208e+46\n", + " 2.616208e+46]\n", + "Rotational Kinetic Energy (J) [3.14564313e+45 3.14564313e+45 3.14564313e+45 3.14564313e+45\n", + " 3.14564313e+45 3.14564313e+45]\n" + ] + } + ], + "source": [ + "#############################################################\n", + "# Print ratios of G.Pot.Energy/B_energy and G.Force/E_force #\n", + "#############################################################\n", + "print(\"Gravitational force to E-force : \", G_EForce)\n", + "print(\"Gravitational energy to B-energy : \",GB_energy_ratio)\n", + "\n", + "\n", + "# Print dimensional values of energies\n", + "print(\"B energy, (J)\",magnetic_energy);\n", + "print(\"Gravitational potential energy, (J)\", GP_energy)\n", + "print(\"Rotational Kinetic Energy (J)\", rotational_KE )\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "application/javascript": [ + "/* Put everything inside the global mpl namespace */\n", + "window.mpl = {};\n", + "\n", + "\n", + "mpl.get_websocket_type = function() {\n", + " if (typeof(WebSocket) !== 'undefined') {\n", + " return WebSocket;\n", + " } else if (typeof(MozWebSocket) !== 'undefined') {\n", + " return MozWebSocket;\n", + " } else {\n", + " alert('Your browser does not have WebSocket support. ' +\n", + " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", + " 'Firefox 4 and 5 are also supported but you ' +\n", + " 'have to enable WebSockets in about:config.');\n", + " };\n", + "}\n", + "\n", + "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", + " this.id = figure_id;\n", + "\n", + " this.ws = websocket;\n", + "\n", + " this.supports_binary = (this.ws.binaryType != undefined);\n", + "\n", + " if (!this.supports_binary) {\n", + " var warnings = document.getElementById(\"mpl-warnings\");\n", + " if (warnings) {\n", + " warnings.style.display = 'block';\n", + " warnings.textContent = (\n", + " \"This browser does not support binary websocket messages. \" +\n", + " \"Performance may be slow.\");\n", + " }\n", + " }\n", + "\n", + " this.imageObj = new Image();\n", + "\n", + " this.context = undefined;\n", + " this.message = undefined;\n", + " this.canvas = undefined;\n", + " this.rubberband_canvas = undefined;\n", + " this.rubberband_context = undefined;\n", + " this.format_dropdown = undefined;\n", + "\n", + " this.image_mode = 'full';\n", + "\n", + " this.root = $('
');\n", + " this._root_extra_style(this.root)\n", + " this.root.attr('style', 'display: inline-block');\n", + "\n", + " $(parent_element).append(this.root);\n", + "\n", + " this._init_header(this);\n", + " this._init_canvas(this);\n", + " this._init_toolbar(this);\n", + "\n", + " var fig = this;\n", + "\n", + " this.waiting = false;\n", + "\n", + " this.ws.onopen = function () {\n", + " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", + " fig.send_message(\"send_image_mode\", {});\n", + " if (mpl.ratio != 1) {\n", + " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", + " }\n", + " fig.send_message(\"refresh\", {});\n", + " }\n", + "\n", + " this.imageObj.onload = function() {\n", + " if (fig.image_mode == 'full') {\n", + " // Full images could contain transparency (where diff images\n", + " // almost always do), so we need to clear the canvas so that\n", + " // there is no ghosting.\n", + " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", + " }\n", + " fig.context.drawImage(fig.imageObj, 0, 0);\n", + " };\n", + "\n", + " this.imageObj.onunload = function() {\n", + " fig.ws.close();\n", + " }\n", + "\n", + " this.ws.onmessage = this._make_on_message_function(this);\n", + "\n", + " this.ondownload = ondownload;\n", + "}\n", + "\n", + "mpl.figure.prototype._init_header = function() {\n", + " var titlebar = $(\n", + " '
');\n", + " var titletext = $(\n", + " '
');\n", + " titlebar.append(titletext)\n", + " this.root.append(titlebar);\n", + " this.header = titletext[0];\n", + "}\n", + "\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._init_canvas = function() {\n", + " var fig = this;\n", + "\n", + " var canvas_div = $('
');\n", + "\n", + " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", + "\n", + " function canvas_keyboard_event(event) {\n", + " return fig.key_event(event, event['data']);\n", + " }\n", + "\n", + " canvas_div.keydown('key_press', canvas_keyboard_event);\n", + " canvas_div.keyup('key_release', canvas_keyboard_event);\n", + " this.canvas_div = canvas_div\n", + " this._canvas_extra_style(canvas_div)\n", + " this.root.append(canvas_div);\n", + "\n", + " var canvas = $('');\n", + " canvas.addClass('mpl-canvas');\n", + " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", + "\n", + " this.canvas = canvas[0];\n", + " this.context = canvas[0].getContext(\"2d\");\n", + "\n", + " var backingStore = this.context.backingStorePixelRatio ||\n", + "\tthis.context.webkitBackingStorePixelRatio ||\n", + "\tthis.context.mozBackingStorePixelRatio ||\n", + "\tthis.context.msBackingStorePixelRatio ||\n", + "\tthis.context.oBackingStorePixelRatio ||\n", + "\tthis.context.backingStorePixelRatio || 1;\n", + "\n", + " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", + "\n", + " var rubberband = $('');\n", + " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", + "\n", + " var pass_mouse_events = true;\n", + "\n", + " canvas_div.resizable({\n", + " start: function(event, ui) {\n", + " pass_mouse_events = false;\n", + " },\n", + " resize: function(event, ui) {\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " stop: function(event, ui) {\n", + " pass_mouse_events = true;\n", + " fig.request_resize(ui.size.width, ui.size.height);\n", + " },\n", + " });\n", + "\n", + " function mouse_event_fn(event) {\n", + " if (pass_mouse_events)\n", + " return fig.mouse_event(event, event['data']);\n", + " }\n", + "\n", + " rubberband.mousedown('button_press', mouse_event_fn);\n", + " rubberband.mouseup('button_release', mouse_event_fn);\n", + " // Throttle sequential mouse events to 1 every 20ms.\n", + " rubberband.mousemove('motion_notify', mouse_event_fn);\n", + "\n", + " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", + " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", + "\n", + " canvas_div.on(\"wheel\", function (event) {\n", + " event = event.originalEvent;\n", + " event['data'] = 'scroll'\n", + " if (event.deltaY < 0) {\n", + " event.step = 1;\n", + " } else {\n", + " event.step = -1;\n", + " }\n", + " mouse_event_fn(event);\n", + " });\n", + "\n", + " canvas_div.append(canvas);\n", + " canvas_div.append(rubberband);\n", + "\n", + " this.rubberband = rubberband;\n", + " this.rubberband_canvas = rubberband[0];\n", + " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", + " this.rubberband_context.strokeStyle = \"#000000\";\n", + "\n", + " this._resize_canvas = function(width, height) {\n", + " // Keep the size of the canvas, canvas container, and rubber band\n", + " // canvas in synch.\n", + " canvas_div.css('width', width)\n", + " canvas_div.css('height', height)\n", + "\n", + " canvas.attr('width', width * mpl.ratio);\n", + " canvas.attr('height', height * mpl.ratio);\n", + " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", + "\n", + " rubberband.attr('width', width);\n", + " rubberband.attr('height', height);\n", + " }\n", + "\n", + " // Set the figure to an initial 600x600px, this will subsequently be updated\n", + " // upon first draw.\n", + " this._resize_canvas(600, 600);\n", + "\n", + " // Disable right mouse context menu.\n", + " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", + " return false;\n", + " });\n", + "\n", + " function set_focus () {\n", + " canvas.focus();\n", + " canvas_div.focus();\n", + " }\n", + "\n", + " window.setTimeout(set_focus, 100);\n", + "}\n", + "\n", + "mpl.figure.prototype._init_toolbar = function() {\n", + " var fig = this;\n", + "\n", + " var nav_element = $('
');\n", + " nav_element.attr('style', 'width: 100%');\n", + " this.root.append(nav_element);\n", + "\n", + " // Define a callback function for later on.\n", + " function toolbar_event(event) {\n", + " return fig.toolbar_button_onclick(event['data']);\n", + " }\n", + " function toolbar_mouse_event(event) {\n", + " return fig.toolbar_button_onmouseover(event['data']);\n", + " }\n", + "\n", + " for(var toolbar_ind in mpl.toolbar_items) {\n", + " var name = mpl.toolbar_items[toolbar_ind][0];\n", + " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", + " var image = mpl.toolbar_items[toolbar_ind][2];\n", + " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", + "\n", + " if (!name) {\n", + " // put a spacer in here.\n", + " continue;\n", + " }\n", + " var button = $('