- Email: [email protected]
Use of a hybrid code for global-scale simulation of the Earth's magnetosphere
Computer Physics Communications 164 (2004) 370–376 www.elsevier.com/locate/cpc
Use of a hybrid code for global-scale simulation of the Earth’s magnetosphere Daniel W. Swift Geophysical Institute, University of Alaska, Fairbanks, AK 99775-7320, USA Available online 21 July 2004
Abstract An understanding of the Earth’s magnetosphere and the magnetospheric substorm requires an ability to simulate kinetic processes in a global context. To this end we have developed a hybrid code with a number of innovative features: (1) The code uses a general curvilinear coordinate grid enabling it to accommodate disparate size scales. (2) The cold, dense plasma inside the plasmapause, where kinetic effects are unimportant, is treated in the fluid approximation. There is a seamless interface between the fluid and kinetic descriptions. (3) The field update is subcycled to the particle update to circumvent the Courant condition with respect to the propagation of Alfvén and whistler mode waves. A two-dimensional version of the code has established a causal link between events in the plasma sheet and the filamentary currents carried by auroral electrons. We are now embarking on the development of a full three-dimensional code for the coupled ionosphere–magnetosphere–magnetosheath system. The simulation domain is comprised of seven discontinuously joined coordinate patches. The Earth’s ionospheric shell constitutes one of the boundary surfaces. 2004 Published by Elsevier B.V. PACS: 94.20.Lr Keywords: Hybrid simulation; Global scale; Magnetosphere
1. Introduction The Earth’s magnetosphere provides the arena for space weather. Dynamical processes in the magnetosphere are responsible for the aurora, geomagnetic storm and substorm phenomena and the electron radiation belt environment. Magnetospheric plasmas exhibit large departures from equilibrium. Global-scale MHD simulations of the magnetosphere have been only partially successful in providing a basic under-
E-mail address: [email protected] (D.W. Swift). 0010-4655/$ – see front matter 2004 Published by Elsevier B.V. doi:10.1016/j.cpc.2004.06.049
standing of space weather. A more complete understanding of the Earth’s magnetosphere and the magnetospheric substorm requires an ability to simulate kinetic processes in a global context. To this end we are in the process of developing a global-scale hybrid code of the coupled ionosphere, magnetosphere, magnetosheath system. The hybrid code provides a full kinetic treatment of the ions, while approximating the electrons as massless fluid. The hybrid code is based on an algorithm was developed by Harned [1]. In a nutshell, the magnetic field in units of ion gyrofrequency is advanced by
D.W. Swift / Computer Physics Communications 164 (2004) 370–376
Faraday’s law ∂B (1) = −∇ × E, ∂t where the electric field in units of ion acceleration is derived from the electron momentum equation E = −ue × B.
(2)
The electron flow velocity is obtained from Ampere’s Law ue = ui −
∇ ×B , αn
α=
4πe2 . mi c 2
than is required to accommodate the Courant condition with respect to the propagation of Alfvén and whistler mode waves in the near-Earth environment. We have typically used 30 time steps to update the fields for every particle time step. An additional modification since the publication of Ref. [2] has been the introduction of the Villasenor–Bunemann [3] exactly conservative current algorithm. Although exact charge conservation is not required in a hybrid code, the algorithm adapts naturally to a structured curvilinear coordinate system.
(3)
The density and ion flow velocity are calculated from the particle moments n= S(x − xk ), nui = vk S(x − xk ). (4) k
371
k
Here S stands for standard PIC weighting, and the summation is over all the particles in the system. The positions and velocities of each particle are advanced by Newton’s laws of motion and the Lorentz force equation. The code used in simulation of the Earth’s magnetosphere has a number of important modifications as described in Swift [2]. (1) The code uses a general curvilinear coordinate grid enabling it to accommodate disparate size scales. This makes it possible to resolve distances of 100 km at the Earth’s ionosphere, while accommodating a simulation domain spanning many tens of Earth radii (RE). (2) The code accommodates both the fully kinetic and MHD description of the ions within the same spatial volume. This provides a seamless interface between the two descriptions. This makes it possible to treat the cold, dense plasma where kinetic effects are unimportant, in the fluid approximation, which provides enormous computational savings, since plasma inside the plasmapause at about 4 RE (Earth radii) constitutes most of the plasma mass. The dense plasma also acts to slow the propagation of Alfvén and whistler mode waves, which mitigates the severe time step constraint in regions where the Earth’s dipole field is strong. (3) The field update is subcycled to the particle update. The most computationally expensive part of the simulation is in the advancement of particle phase space positions. Since most of the kinetic particles are beyond 4 RE, they can be advanced with a longer time step
2. Results from two-dimensional simulation This section will show the results of a simulation in the midnight meridian plane of the Earth’s magnetosphere. The simulation domain extends from the polar axis of the Earth to 36 Earth radii in the antisunward direction. The usual GSM coordinates are used with x being in the sunward direction, y in the dawn-to-dusk direction and z aligned along the geomagnetic polar axes. Distances are measured in Earth radii. The run is initiated with an extended magnetotail in which the sunward-directed Maxwell stresses on the plasma sheet exceeded the antisunward pressure gradient forces. Fig. 1 shows the results of the simulation at two different times. The top panels show the field line traces. The unbalanced Maxwell stress across the plasma sheet accelerates plasma earthward, as shown in the middle panels. The flow velocities reach 500 km/s, about two-thirds the Alfvén rate. The flow is stopped and diverted by the Earth’s dipole field. The bottom panel shows the dawn-dusk component of the magnetic field. The fluctuation amplitudes are about 12 nT, just somewhat less than the 15 nT magnitude of the asymptotic lobe field. It can be seen that turbulence is generated in region where the flow is stopped. More detailed diagnostics shown by Swift and Lin [4] show multiple peak ion velocity distributions caused by the interpenetration of the earthward flowing and the ambient plasma. The run was carried out to t = 600 s, and as the run progressed the turbulent region retreated antisunward. Fig. 2(a), shows a more detailed view of the dawndusk component of the magnetic field. An animation shows that the turbulent region radiates shear Alfvén waves that propagate along magnetic field lines toward
372
D.W. Swift / Computer Physics Communications 164 (2004) 370–376
Fig. 1. The panels show meridian plane views. The small semicircle on the right hand side of the panels represents the region occupied by the Earth. The Sun is to the right, and distances are measured in Earth radii. The top panels show field line traces at the times indicated in seconds. The middle panels show plasma flow. The bottom panels show the dawn-dusk component of the magnetic field, which is normal to the plane of simulation.
the high latitude ionosphere. Panel (b) shows a close up view of the field-aligned component of the current, obtained by differentiating the y-component of
the magnetic field. These filamentary currents of the type associated with proton and auroral electron precipitation are propagated by the Alfvén waves. The an-
D.W. Swift / Computer Physics Communications 164 (2004) 370–376
Fig. 2. A more detailed view of the dawn-dusk component of the magnetic field is shown in panel (a). Panel (b) shows a close-up view of the field-aligned currents in the plane of the simulation in the near-Earth region. The current is derived from the curl of the dawn-to-dusk, or y-component of the magnetic field.
tisunward retreat of the region of flow breaking and plasma sheet turbulence causes a poleward progression of the filamentary currents at the Earth’s surface. This is consistent with the poleward progression of auroral arcs during a substorm expansion. Moreover, recent satellite observations of Wygant et al. [5] in the midway region between the plasma sheet and auroral ionosphere have shown large amplitude Alfvén waves carrying Poynting fluxes considerably in excess of that necessary to account for auroral fluxes precipitating into the atmosphere. The two-dimensional hybrid code has for the first time shown a causal link between fast flows in the plasma sheet [6] caused by the relaxation of an extended magnetotail and the field-aligned currents above the auroral ionosphere. The generation of the turbulence shown in Figs. 1 and 2 is a kinetic process. The ability to simulate the generation and propagation of the Alfvén waves requires a global setting.
373
The two-dimensional simulation has a number of limitations. It is unable to capture the westward electrojet associated with auroral substorms, since it involves structure in the dawn-dusk direction. The two-dimensional simulations show a pileup of hot plasma in the stagnation region where the earthward plasma flow is broken. In three dimensions this plasma cloud flows around to the afternoon side of the magnetosphere. In three dimensions, the cloud can also polarize to allow it to penetrate into the inner magnetosphere to form the ring current, which is responsible for the main phase, storm-time magnetic field decrease at low latitudes. However, since the two-dimensional simulation has shown the Alfvén wave emitting turbulence to be generated on the outer boundary of the plasma-pileup region, we would expect the tree-dimensional code would also show the generation of turbulence where the earthward fast flow is broken. Also, with the two-dimensional code, we are unable to simulate the stretching of the tail field lines during the period, known as the growth phase, leading up to the substorm expansion. We believe the stretching process involves the ability of plasmoids associated with bursty bulk flow events to generate dawn-dusk electrostatic fields. The successes and limitations of the two-dimensional model provide the motivation to undertake the much larger task of developing a global-scale three-dimensional that includes ion kinetics.
3. The three-dimensional model The three-dimensional model is built around seven discontinuously joined coordinate patches. Multiple coordinate patches are used in order to bring the coordinate grid to the surface of the Earth while avoiding coordinate singularities. Fig. 3(a) shows a meridian cut of the grid. The sun is to the right. The outer boundaries at the sunward, poleward and dawn and dusk patches are at about 15 RE. The rectangular tail patch extends to 40 RE in the antisunward direction. The inner boundary is the Earth’s ionosphere at 1 RE, where the grid resolution is less than 100 km. The six coordinate patches bounding on the Earth are shaped like a pyramid with a spherical cavity at the apex, as shown in panel (b) of Fig. 3. Two of the coordinate patches not shown are in the dawn and dusk sec-
374
D.W. Swift / Computer Physics Communications 164 (2004) 370–376
Fig. 3. A meridian cut of the three-dimensional grid is shown in panel (a). The grid extends to 15 Earth radii on the sunward, north, south, dawn and dusk sides and extends beyond 40 Earth radii on the night side. The Earth is shown as the circular blank spot. Panel (b) shows a perspective view of one of the coordinate patches that bounds on the surface of the Earth.
tors. Notice the sunward patch has a greater density of grid points beyond about 9 RE to resolve the bow shock and magnetopause. Also, there is a greater density of grid points near the equatorial regions of the antisunward and tail patches to resolve important features of the plasma sheet. No attempt has been made
to match grid lines across coordinate patches; however, great care has been taken to ensure continuity of the density of coordinate surfaces across patch boundaries to avoid impedance mismatch for waves propagating across patch boundaries. The coordinates are specified by a table giving the x, y, z-locations of the
D.W. Swift / Computer Physics Communications 164 (2004) 370–376
grid points. All the geometric coefficients necessary for field operations and moving particles in a curvilinear coordinate system are derived from the table. The grid can be easily modified without having to do any reprogramming by simply changing the underlying table specifying the grid point locations. The grids are set up so there is a half grid width overlap between coordinate patches. Field data from one coordinate patch is interpolated, in both position and direction, onto the overlapping grid points of the neighboring patch. The divergence of the electric field, E, and the magnetic field, B, are on nested grids, with components of the fields residing on the cell faces. Boundary conditions are imposed on the tangential components of the electric field and the normal component of the magnetic field. These components must be interpolated from neighboring patches. Similarly, the particle curvilinear position has to be interpolated from one coordinate patch to another as the particle crosses patch boundaries. No conversion of the particle velocity is required, as the velocity is kept in Cartesian coordinates. The field matching procedures have been tested using test fields with constant curl. The test fields are generated in Cartesian coordinates, then converted to curvilinear, and the curl is taken in curvilinear coordinates. The curl is converted back to Cartesian and compared with the constant Cartesian value. The error is less than 0.1% in the interior of coordinate patches and less than 1% at patch boundaries and edges. Experience with the two-dimensional simulations has shown this accuracy to be sufficient. The particle mover has been tested by tracing the trajectory of a force-free particle across coordinate patch boundaries. Trajectory plots have shown no discernable deviation from straightline trajectories. The code is being set up for a scalable, distributed memory parallel computer. In order to have some flexibility in the assignment of processors, we have taken the additional step of making it possible to independently subdivide each coordinate patch into a number of different segments to accommodate a greater number of processors. The grid points and particles within each segment are assigned to a single processor. There is no requirement for matching segment boundaries across patch boundaries. The configuration of segments within each processor is specified by a set of indices giving the beginning and ending coordinate
375
indices of each segment within each patch. The segment configuration within each coordinate patch can be altered without reprogramming by changing the indices defining its extent. Passing of field and particle data within each coordinate patch occurs directly between the segments. Segments having a face on a coordinate patch boundary pass its guard cell data to a processor designated to assemble data for an entire patch face. Individual segments from the adjoining coordinate patch then access their individual pieces of surface data from the processor that had assembled the data. This two-step process provides the flexibility that enables the coordinate patches to be independently subdivided. External boundary conditions are handled in a similar way. A single processor develops the boundary data for an external face of an entire coordinate patch, and the individual segments access the boundary data for their external boundary face. An important “external boundary” is that presented by the partially conducting ionosphere. The normal component of the magnetic field is assumed to be that of the Earth’s dipole. The tangential component of the electric field is derived from the condition that the sum of the normal component of the current and the horizontal, height-integrated current flowing in the ionospheric shell be divergenceless. The normal component of the current is derived from the curl of the tangential components of the magnetic field a half-grid layer above the ionospheric shell. The horizontal component of the electric field is assumed to be described by a potential. The conductivity coefficients are derived from the collisional coupling between the ions and electrons, and the neutral atmosphere. Co-rotation in the frame fixed with respect to the Earth–sun line is imposed by having the neutral atmosphere rotate with the Earth. The ions and electrons within the ionosphere are driven by collisions with the far more abundant neutrals. The conductivity dyadic is anisotropic because of the Hall term. The derived conductivity coefficients have large spatial variation depending on the ionization density in the conducting layer. The potential is calculated on the 102 × 102 polar grids of the north and south polar patches. The electric fields on the lower latitude grids will be derived from the co-rotation potential. The entire system is driven by the solar wind and the interplanetary magnetic field (IMF) carried by
376
D.W. Swift / Computer Physics Communications 164 (2004) 370–376
the solar wind. Responses of the magnetosphere to changes in the IMF are well known. The code will be tested by subjecting the magnetosphere to changes in IMF conditions to see how well the model responses align with observations.
rora and magnetospheric dynamics. Further progress requires three-dimensional simulations of the coupled ionosphere, magnetosphere and solar wind. We are now in the process of building a global scale model consisting of seven discontinuously joined coordinate patches. This model is being designed to run on a scalable massively parallel computer.
4. Summary Use of a two-dimensional hybrid code in a globalscale simulation of the Earth’s magnetosphere has revealed new physics. Namely, the code has shown a causal connection between the relaxation of an extended magnetotail field, fast earthward flows, generation of large amplitude Alfvén waves and production of aurora. Simulations involving ion kinetics in a global context is essential to understanding the au-
References [1] D.S. Harned, J. Comput. Phys. 47 (1982) 452. [2] D.W. Swift, J. Comput. Phys. 126 (1996) 109. [3] J. Villasenor, O. Bunemann, Comput. Phys. Commun. 2 (1982) 306. [4] D.W. Swift, Y. Lin, JASTP 63 (2001) 683. [5] J.R. Wygant, et al., J. Geophys. Res. 105 (2000) 18,765. [6] J.-H. Shue, et al., J. Geophys. Res. 108 (2003), 10.1029/ 2002JA009794.
Use of a hybrid code for global-scale simulation of the Earth’s magnetosphere Daniel W. Swift Geophysical Institute, University of Alaska, Fairbanks, AK 99775-7320, USA Available online 21 July 2004
Abstract An understanding of the Earth’s magnetosphere and the magnetospheric substorm requires an ability to simulate kinetic processes in a global context. To this end we have developed a hybrid code with a number of innovative features: (1) The code uses a general curvilinear coordinate grid enabling it to accommodate disparate size scales. (2) The cold, dense plasma inside the plasmapause, where kinetic effects are unimportant, is treated in the fluid approximation. There is a seamless interface between the fluid and kinetic descriptions. (3) The field update is subcycled to the particle update to circumvent the Courant condition with respect to the propagation of Alfvén and whistler mode waves. A two-dimensional version of the code has established a causal link between events in the plasma sheet and the filamentary currents carried by auroral electrons. We are now embarking on the development of a full three-dimensional code for the coupled ionosphere–magnetosphere–magnetosheath system. The simulation domain is comprised of seven discontinuously joined coordinate patches. The Earth’s ionospheric shell constitutes one of the boundary surfaces. 2004 Published by Elsevier B.V. PACS: 94.20.Lr Keywords: Hybrid simulation; Global scale; Magnetosphere
1. Introduction The Earth’s magnetosphere provides the arena for space weather. Dynamical processes in the magnetosphere are responsible for the aurora, geomagnetic storm and substorm phenomena and the electron radiation belt environment. Magnetospheric plasmas exhibit large departures from equilibrium. Global-scale MHD simulations of the magnetosphere have been only partially successful in providing a basic under-
E-mail address: [email protected] (D.W. Swift). 0010-4655/$ – see front matter 2004 Published by Elsevier B.V. doi:10.1016/j.cpc.2004.06.049
standing of space weather. A more complete understanding of the Earth’s magnetosphere and the magnetospheric substorm requires an ability to simulate kinetic processes in a global context. To this end we are in the process of developing a global-scale hybrid code of the coupled ionosphere, magnetosphere, magnetosheath system. The hybrid code provides a full kinetic treatment of the ions, while approximating the electrons as massless fluid. The hybrid code is based on an algorithm was developed by Harned [1]. In a nutshell, the magnetic field in units of ion gyrofrequency is advanced by
D.W. Swift / Computer Physics Communications 164 (2004) 370–376
Faraday’s law ∂B (1) = −∇ × E, ∂t where the electric field in units of ion acceleration is derived from the electron momentum equation E = −ue × B.
(2)
The electron flow velocity is obtained from Ampere’s Law ue = ui −
∇ ×B , αn
α=
4πe2 . mi c 2
than is required to accommodate the Courant condition with respect to the propagation of Alfvén and whistler mode waves in the near-Earth environment. We have typically used 30 time steps to update the fields for every particle time step. An additional modification since the publication of Ref. [2] has been the introduction of the Villasenor–Bunemann [3] exactly conservative current algorithm. Although exact charge conservation is not required in a hybrid code, the algorithm adapts naturally to a structured curvilinear coordinate system.
(3)
The density and ion flow velocity are calculated from the particle moments n= S(x − xk ), nui = vk S(x − xk ). (4) k
371
k
Here S stands for standard PIC weighting, and the summation is over all the particles in the system. The positions and velocities of each particle are advanced by Newton’s laws of motion and the Lorentz force equation. The code used in simulation of the Earth’s magnetosphere has a number of important modifications as described in Swift [2]. (1) The code uses a general curvilinear coordinate grid enabling it to accommodate disparate size scales. This makes it possible to resolve distances of 100 km at the Earth’s ionosphere, while accommodating a simulation domain spanning many tens of Earth radii (RE). (2) The code accommodates both the fully kinetic and MHD description of the ions within the same spatial volume. This provides a seamless interface between the two descriptions. This makes it possible to treat the cold, dense plasma where kinetic effects are unimportant, in the fluid approximation, which provides enormous computational savings, since plasma inside the plasmapause at about 4 RE (Earth radii) constitutes most of the plasma mass. The dense plasma also acts to slow the propagation of Alfvén and whistler mode waves, which mitigates the severe time step constraint in regions where the Earth’s dipole field is strong. (3) The field update is subcycled to the particle update. The most computationally expensive part of the simulation is in the advancement of particle phase space positions. Since most of the kinetic particles are beyond 4 RE, they can be advanced with a longer time step
2. Results from two-dimensional simulation This section will show the results of a simulation in the midnight meridian plane of the Earth’s magnetosphere. The simulation domain extends from the polar axis of the Earth to 36 Earth radii in the antisunward direction. The usual GSM coordinates are used with x being in the sunward direction, y in the dawn-to-dusk direction and z aligned along the geomagnetic polar axes. Distances are measured in Earth radii. The run is initiated with an extended magnetotail in which the sunward-directed Maxwell stresses on the plasma sheet exceeded the antisunward pressure gradient forces. Fig. 1 shows the results of the simulation at two different times. The top panels show the field line traces. The unbalanced Maxwell stress across the plasma sheet accelerates plasma earthward, as shown in the middle panels. The flow velocities reach 500 km/s, about two-thirds the Alfvén rate. The flow is stopped and diverted by the Earth’s dipole field. The bottom panel shows the dawn-dusk component of the magnetic field. The fluctuation amplitudes are about 12 nT, just somewhat less than the 15 nT magnitude of the asymptotic lobe field. It can be seen that turbulence is generated in region where the flow is stopped. More detailed diagnostics shown by Swift and Lin [4] show multiple peak ion velocity distributions caused by the interpenetration of the earthward flowing and the ambient plasma. The run was carried out to t = 600 s, and as the run progressed the turbulent region retreated antisunward. Fig. 2(a), shows a more detailed view of the dawndusk component of the magnetic field. An animation shows that the turbulent region radiates shear Alfvén waves that propagate along magnetic field lines toward
372
D.W. Swift / Computer Physics Communications 164 (2004) 370–376
Fig. 1. The panels show meridian plane views. The small semicircle on the right hand side of the panels represents the region occupied by the Earth. The Sun is to the right, and distances are measured in Earth radii. The top panels show field line traces at the times indicated in seconds. The middle panels show plasma flow. The bottom panels show the dawn-dusk component of the magnetic field, which is normal to the plane of simulation.
the high latitude ionosphere. Panel (b) shows a close up view of the field-aligned component of the current, obtained by differentiating the y-component of
the magnetic field. These filamentary currents of the type associated with proton and auroral electron precipitation are propagated by the Alfvén waves. The an-
D.W. Swift / Computer Physics Communications 164 (2004) 370–376
Fig. 2. A more detailed view of the dawn-dusk component of the magnetic field is shown in panel (a). Panel (b) shows a close-up view of the field-aligned currents in the plane of the simulation in the near-Earth region. The current is derived from the curl of the dawn-to-dusk, or y-component of the magnetic field.
tisunward retreat of the region of flow breaking and plasma sheet turbulence causes a poleward progression of the filamentary currents at the Earth’s surface. This is consistent with the poleward progression of auroral arcs during a substorm expansion. Moreover, recent satellite observations of Wygant et al. [5] in the midway region between the plasma sheet and auroral ionosphere have shown large amplitude Alfvén waves carrying Poynting fluxes considerably in excess of that necessary to account for auroral fluxes precipitating into the atmosphere. The two-dimensional hybrid code has for the first time shown a causal link between fast flows in the plasma sheet [6] caused by the relaxation of an extended magnetotail and the field-aligned currents above the auroral ionosphere. The generation of the turbulence shown in Figs. 1 and 2 is a kinetic process. The ability to simulate the generation and propagation of the Alfvén waves requires a global setting.
373
The two-dimensional simulation has a number of limitations. It is unable to capture the westward electrojet associated with auroral substorms, since it involves structure in the dawn-dusk direction. The two-dimensional simulations show a pileup of hot plasma in the stagnation region where the earthward plasma flow is broken. In three dimensions this plasma cloud flows around to the afternoon side of the magnetosphere. In three dimensions, the cloud can also polarize to allow it to penetrate into the inner magnetosphere to form the ring current, which is responsible for the main phase, storm-time magnetic field decrease at low latitudes. However, since the two-dimensional simulation has shown the Alfvén wave emitting turbulence to be generated on the outer boundary of the plasma-pileup region, we would expect the tree-dimensional code would also show the generation of turbulence where the earthward fast flow is broken. Also, with the two-dimensional code, we are unable to simulate the stretching of the tail field lines during the period, known as the growth phase, leading up to the substorm expansion. We believe the stretching process involves the ability of plasmoids associated with bursty bulk flow events to generate dawn-dusk electrostatic fields. The successes and limitations of the two-dimensional model provide the motivation to undertake the much larger task of developing a global-scale three-dimensional that includes ion kinetics.
3. The three-dimensional model The three-dimensional model is built around seven discontinuously joined coordinate patches. Multiple coordinate patches are used in order to bring the coordinate grid to the surface of the Earth while avoiding coordinate singularities. Fig. 3(a) shows a meridian cut of the grid. The sun is to the right. The outer boundaries at the sunward, poleward and dawn and dusk patches are at about 15 RE. The rectangular tail patch extends to 40 RE in the antisunward direction. The inner boundary is the Earth’s ionosphere at 1 RE, where the grid resolution is less than 100 km. The six coordinate patches bounding on the Earth are shaped like a pyramid with a spherical cavity at the apex, as shown in panel (b) of Fig. 3. Two of the coordinate patches not shown are in the dawn and dusk sec-
374
D.W. Swift / Computer Physics Communications 164 (2004) 370–376
Fig. 3. A meridian cut of the three-dimensional grid is shown in panel (a). The grid extends to 15 Earth radii on the sunward, north, south, dawn and dusk sides and extends beyond 40 Earth radii on the night side. The Earth is shown as the circular blank spot. Panel (b) shows a perspective view of one of the coordinate patches that bounds on the surface of the Earth.
tors. Notice the sunward patch has a greater density of grid points beyond about 9 RE to resolve the bow shock and magnetopause. Also, there is a greater density of grid points near the equatorial regions of the antisunward and tail patches to resolve important features of the plasma sheet. No attempt has been made
to match grid lines across coordinate patches; however, great care has been taken to ensure continuity of the density of coordinate surfaces across patch boundaries to avoid impedance mismatch for waves propagating across patch boundaries. The coordinates are specified by a table giving the x, y, z-locations of the
D.W. Swift / Computer Physics Communications 164 (2004) 370–376
grid points. All the geometric coefficients necessary for field operations and moving particles in a curvilinear coordinate system are derived from the table. The grid can be easily modified without having to do any reprogramming by simply changing the underlying table specifying the grid point locations. The grids are set up so there is a half grid width overlap between coordinate patches. Field data from one coordinate patch is interpolated, in both position and direction, onto the overlapping grid points of the neighboring patch. The divergence of the electric field, E, and the magnetic field, B, are on nested grids, with components of the fields residing on the cell faces. Boundary conditions are imposed on the tangential components of the electric field and the normal component of the magnetic field. These components must be interpolated from neighboring patches. Similarly, the particle curvilinear position has to be interpolated from one coordinate patch to another as the particle crosses patch boundaries. No conversion of the particle velocity is required, as the velocity is kept in Cartesian coordinates. The field matching procedures have been tested using test fields with constant curl. The test fields are generated in Cartesian coordinates, then converted to curvilinear, and the curl is taken in curvilinear coordinates. The curl is converted back to Cartesian and compared with the constant Cartesian value. The error is less than 0.1% in the interior of coordinate patches and less than 1% at patch boundaries and edges. Experience with the two-dimensional simulations has shown this accuracy to be sufficient. The particle mover has been tested by tracing the trajectory of a force-free particle across coordinate patch boundaries. Trajectory plots have shown no discernable deviation from straightline trajectories. The code is being set up for a scalable, distributed memory parallel computer. In order to have some flexibility in the assignment of processors, we have taken the additional step of making it possible to independently subdivide each coordinate patch into a number of different segments to accommodate a greater number of processors. The grid points and particles within each segment are assigned to a single processor. There is no requirement for matching segment boundaries across patch boundaries. The configuration of segments within each processor is specified by a set of indices giving the beginning and ending coordinate
375
indices of each segment within each patch. The segment configuration within each coordinate patch can be altered without reprogramming by changing the indices defining its extent. Passing of field and particle data within each coordinate patch occurs directly between the segments. Segments having a face on a coordinate patch boundary pass its guard cell data to a processor designated to assemble data for an entire patch face. Individual segments from the adjoining coordinate patch then access their individual pieces of surface data from the processor that had assembled the data. This two-step process provides the flexibility that enables the coordinate patches to be independently subdivided. External boundary conditions are handled in a similar way. A single processor develops the boundary data for an external face of an entire coordinate patch, and the individual segments access the boundary data for their external boundary face. An important “external boundary” is that presented by the partially conducting ionosphere. The normal component of the magnetic field is assumed to be that of the Earth’s dipole. The tangential component of the electric field is derived from the condition that the sum of the normal component of the current and the horizontal, height-integrated current flowing in the ionospheric shell be divergenceless. The normal component of the current is derived from the curl of the tangential components of the magnetic field a half-grid layer above the ionospheric shell. The horizontal component of the electric field is assumed to be described by a potential. The conductivity coefficients are derived from the collisional coupling between the ions and electrons, and the neutral atmosphere. Co-rotation in the frame fixed with respect to the Earth–sun line is imposed by having the neutral atmosphere rotate with the Earth. The ions and electrons within the ionosphere are driven by collisions with the far more abundant neutrals. The conductivity dyadic is anisotropic because of the Hall term. The derived conductivity coefficients have large spatial variation depending on the ionization density in the conducting layer. The potential is calculated on the 102 × 102 polar grids of the north and south polar patches. The electric fields on the lower latitude grids will be derived from the co-rotation potential. The entire system is driven by the solar wind and the interplanetary magnetic field (IMF) carried by
376
D.W. Swift / Computer Physics Communications 164 (2004) 370–376
the solar wind. Responses of the magnetosphere to changes in the IMF are well known. The code will be tested by subjecting the magnetosphere to changes in IMF conditions to see how well the model responses align with observations.
rora and magnetospheric dynamics. Further progress requires three-dimensional simulations of the coupled ionosphere, magnetosphere and solar wind. We are now in the process of building a global scale model consisting of seven discontinuously joined coordinate patches. This model is being designed to run on a scalable massively parallel computer.
4. Summary Use of a two-dimensional hybrid code in a globalscale simulation of the Earth’s magnetosphere has revealed new physics. Namely, the code has shown a causal connection between the relaxation of an extended magnetotail field, fast earthward flows, generation of large amplitude Alfvén waves and production of aurora. Simulations involving ion kinetics in a global context is essential to understanding the au-
References [1] D.S. Harned, J. Comput. Phys. 47 (1982) 452. [2] D.W. Swift, J. Comput. Phys. 126 (1996) 109. [3] J. Villasenor, O. Bunemann, Comput. Phys. Commun. 2 (1982) 306. [4] D.W. Swift, Y. Lin, JASTP 63 (2001) 683. [5] J.R. Wygant, et al., J. Geophys. Res. 105 (2000) 18,765. [6] J.-H. Shue, et al., J. Geophys. Res. 108 (2003), 10.1029/ 2002JA009794.