{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Mapping and Integrating Along Field Lines\n\nThis example demonstrates how to use the :func:`~mapflpy.scripts.map_pt_forward` and\n:func:`~mapflpy.scripts.map_pt_backward` functions to efficiently map large numbers\nof field lines on a grid from one spherical surface to another. These are based off\na very general method for mapping field lines in parallel:\n:func:`~mapflpy.scripts.map_field_lines_in_parallel`.\n\nThis example also illustrates a related task: integrating scalar quantities along field lines.\nThis can be used to compute scalar averages or other quantities. Weighting by the flux-tube\narea expansion is also possible. If no scalars are provided, the field-line length\nis provided by default.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import os\nimport tempfile\n\nimport numpy as np\nfrom psi_io import wrhdf_3d, interpolate_positions_from_hdf\nfrom psi_data import fetch_mas_data\n\n\nimport matplotlib\nimport matplotlib.pyplot as plt\n\nfrom mapflpy.scripts import map_pt_forward, map_pt_backward" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Load in the magnetic field and scalar field files used for this example.\nThey are from a CORHEL-MAS thermodynamic MHD calculation for CR2282.\nBecuase they are not part of the standard datasets in mapflpy or psi-io\n(yet), we fetch them manually and place them in the default cache location.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "files = fetch_mas_data(domains=\"cor\", variables=\"br,bt,bp,t\")\nmagnetic_field_files = files.cor_br, files.cor_bt, files.cor_bp\nscalar_field_file = files.cor_t" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we compute a basic map by mapping field lines forward into the\nvolume from a 2D grid on the surface and recording their endpoint locations.\nThe mapping function accepts 1D arrays that specify the phi (longitude) and\ntheta (co-latitude) grid in radians. It accepts the same mapfl options and keywords\nas the tracing routines (e.g. :class:`~mapflpy.tracer.Tracer`\nor :func:`~mapflpy.scripts.run_forward_tracing`). Here the grid resolution is\n1 degree, which should map relatively quickly on four processors (5-20 seconds).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# define the grid\nradius = 1.0\nt = np.linspace(0, np.pi, 181)\np = np.linspace(0, 2*np.pi, 361)\n\n# compute the mapping\nnproc = 4\nmapping = map_pt_forward(*magnetic_field_files, p, t, radius=radius, nproc=nproc)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we can inspect the :class:`~mapflpy.globals.Mapping` object, which\ncontains the endpoint `r`, `t`, `p` values (2D arrays on the mapping grid in this case).\nHere we plot the endpoint radius, which indicates if a field line is open (30.0)\nor closed (1.0).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ax = plt.figure().add_subplot()\nax.pcolormesh(np.rad2deg(p), 90 - np.rad2deg(t), mapping.r,\n cmap='gray', shading='gouraud', clim=(1, 30))\nax.set_aspect(\"equal\", adjustable=\"box\")\nax.set_title('Map Forward: Endpoint Radius')\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The mapping object also contains the field line length as the default\n`integral` field in units of solar radii. Here we plot it in log space\nto illustrate the dynamic range\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ax = plt.figure().add_subplot()\nax.pcolormesh(np.rad2deg(p), 90 - np.rad2deg(t), np.log10(mapping.integral),\n cmap='gist_ncar', shading='gouraud', clim=(-2, 2))\nax.set_aspect(\"equal\", adjustable=\"box\")\nax.set_title('Map Forward: Field Line Length')\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Instead of length we can integrate a 3D scalar along the field\nby setting mapfl accordingly. Currently, the scalar field must\nbe passed as a path to a 3D HDF5 file. In this example we integrate\ntemperature with a second mapping call. The average temperature\nalong the magnetic field can be determined by dividing this integral\nby the length integral.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mapping2 = map_pt_forward(*magnetic_field_files, p, t, radius=radius, nproc=nproc,\n integrate_along_fl_=True, scalar_input_file_=scalar_field_file)\n\n# compute the average temperature in code units by dividing the two integrals\navg_t_mas = mapping2.integral/mapping.integral\n\n# convert it to MK using the MAS normalizations\nfrom psi_io.units import FN_T\nimport astropy.units as u\n\navg_t_mk = (avg_t_mas*FN_T).to(u.MK)\nax = plt.figure().add_subplot()\nax.pcolormesh(np.rad2deg(p), 90 - np.rad2deg(t), avg_t_mk.value,\n cmap='rainbow', shading='gouraud', clim=(0.5, 2.5))\nax.set_aspect(\"equal\", adjustable=\"box\")\nax.set_title('Map Forward: Average Field Line Temperature')\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also weight the integrals by the relative flux-tube area\nby adding 1/B weighting. This allows one to construct other types of\ninteresting scalar integrals. Here we use it to compute the relative\nfield line volume: $V_{FL} = \\int_0^L A ds = A_0 B_0 \\int_0^L \\frac{1}{B}ds = |d\\vec{A}\\cdot \\vec{B}_0| \\int_0^L \\frac{1}{B}ds= r^2d\\Omega |B_r| \\int_0^L \\frac{1}{B}ds$\nWe do this by turning on `weight_integral_by_area_` and passing\na scalar field file that is equal to 1 at all points.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# create an interior p, t mesh (no poles) to make the area calculation simpler\np_i = 0.5*(p[1:] + p[0:-1])\nt_i = 0.5*(t[1:] + t[0:-1])\n\n# map foward from 1.0 to 2.0\ndomain_r_min = 1.0\ndomain_r_max = 2.0\n\n# Build an all-ones scalar field for the 1/B area weighting. The value is\n# constant, so a coarse (r, theta, phi) grid suffices \u2014 only the scales (which\n# must bracket the traced volume) and the 3D-with-scales format matter. mapfl\n# reads the scalar from disk, so we write it to a temporary file.\nones_r = np.array([domain_r_min, domain_r_max])\nones_t = np.array([0.0, np.pi])\nones_p = np.array([0.0, 2*np.pi])\nones_field = np.ones((ones_p.size, ones_t.size, ones_r.size)) # C-order: (nphi, ntheta, nr)\ndummy_file = os.path.join(tempfile.gettempdir(), 'mapflpy_dummy_ones_3d.h5')\nwrhdf_3d(dummy_file, ones_r, ones_t, ones_p, ones_field)\nmappping_area_fwd = map_pt_forward(*magnetic_field_files, p_i, t_i, radius=domain_r_min, nproc=nproc,\n domain_r_min_=domain_r_min, domain_r_max_=domain_r_max,\n integrate_along_fl_=True, weight_integral_by_area_=True,\n scalar_input_file_=dummy_file)\n\n# compute the area of each gridcell, adjusting the area factor at the pole for the half cell\n# (use np.gradient with care, this is a simple example where it is fine).\ndp = np.gradient(p_i)\ndt = np.gradient(t_i)\ncell_omega = np.einsum('i,j->ij', dt*np.sin(t_i), dp)\n\n# get br at this surface by interpolating to the 2D grid of r, t, p positions\np2d, t2d = np.meshgrid(p_i, t_i)\nones2d = np.ones_like(p2d)\nbr_lower = interpolate_positions_from_hdf(files.cor_br, ones2d*domain_r_min, t2d, p2d)\n\n# Compute the volume\nvol_fwd = domain_r_min**2*cell_omega*np.abs(br_lower)*mappping_area_fwd.integral\n\nax = plt.figure().add_subplot()\nax.pcolormesh(np.rad2deg(p_i), 90 - np.rad2deg(t_i), np.log10(np.clip(vol_fwd, min=1e-5)),\n cmap='terrain', shading='gouraud', clim=(-5, -2))\nax.set_aspect(\"equal\", adjustable=\"box\")\nax.set_title('Map Forward: Field Line Volume [$R_S^3$]')\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lastly, we can confirm (for fun) that the volume calculation is correct by\ncomparing the volume of the fowards and backwards maps to the analytic volume\nof a sphere.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# map backwards from radius=rmax\nmapping_area_bwd = map_pt_backward(*magnetic_field_files, p_i, t_i, radius=domain_r_max, nproc=nproc,\n domain_r_min_=domain_r_min, domain_r_max_=domain_r_max,\n integrate_along_fl_=True, weight_integral_by_area_=True,\n scalar_input_file_=dummy_file)\n\nbr_upper = interpolate_positions_from_hdf(files.cor_br, ones2d*domain_r_max, t2d, p2d)\nvol_bwd = domain_r_max**2*cell_omega*np.abs(br_upper)*mapping_area_bwd.integral\n\n# the total volume should be the average of the two mapping's total volumes\n# (all closed, open, and disconnected field lines only counted once)\nvol_avg = 0.5*(np.sum(vol_fwd) + np.sum(vol_bwd))\n\n# compare to the analytic volume\n# NOTE: This agreement improves with more points in the mapping.\nvol_analytic = 4./3.*np.pi*(domain_r_max**3 - domain_r_min**3)\nprint(f'### Volume check for {len(p_i)}x{len(t_i)} mappings from r = {domain_r_min} - {domain_r_max} Rs')\nprint(f' Analytic Volume: {vol_analytic:.4f}')\nprint(f' Avg Total FL Volume: {vol_avg:.4f}')\nprint(f' Percentage error: {(vol_analytic - vol_avg)/vol_analytic*100:.2f} %')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.13.13" } }, "nbformat": 4, "nbformat_minor": 0 }