User_Manual
This User's Manual is intended to help scientists and programmers use and develop the FLUX source code. It is itself under development. As of Spring 2008, it is primarily a description and design document outlining the major functionality of the code. Sections on how to use the code are in progress.
FLUX currently calculates only the Lorenz force, which makes it a nonlinear force-free field solver (NLFFF solver) with fixed topology. This inverts the boundary/initial condition of most existing codes: most NLFFF codes accept a field angle at the boundaries and return a topology. FLUX instead accepts a topology and returns a field angle at the boundary.
This makes it useful for studying energy release in flux systems with a prescribed topology, or where the topology (but not necessarily field angle) has been matched to some other set of observations, e.g. EUV coronal images.
"How do you solve problems that you don't already know the answer to?" - a solar MHD modeler, 2005
FLUX is a bit unusual in that it is adapted for a different set of boundary conditions than most other MHD or NLFFF codes. The input boundary condition is a topology rather than a field angle at the boundary. Codes with fixed or arbitrary grids model field angle and strength explicitly everywhere; to visualize the topology, one shoots field lines through the calculated fields in the simulation arena. Field direction and strength at arbitrary locations are secondary variables for FLUX, and must be interpolated from the surrounding fluxon geometry. Field angle at the boundary is an output (dependent) variable of the code, rather than a specified initial condition.
Because FLUX conserves topology explicitly, in non-reconnective simulations you must specify exactly what initial topology you want for the simulation. The main use of FLUX as a standalone model is tracking ideal evolution and/or injecting reconnection in a controlled manner. As a semi-empirical model, FLUX can be driven by feature-tracked or extracted velocitiy maps at the surface, but requires an initial topology. FLUX can track gradual buildup of stress in tangled or non-trivial systems that are not essily accessible to more conventional codes.
Several tools exist to input initial topologies into FLUX. The principal tools used so far are:
You can of course hand-tweak fluxon geometries a priori. The simplest ways to to do this are to generate a fragment of an ASCII-format FLUX file with printf or the equivalent, or (in Perl/PDL) to generate a Perl list of 3×ni PDLs that each describe the vertices in one fluxon. There is an autoload routine, make_world(), that will convert the list of PDLs to an ASCII FLUX World string. That string can, in turn, be made into a Flux::World object with str2world.
You can define a collection of a priori domains by declaring a collection of semicircular loops, each of which has simple connectivity. This is the method used by Rachmeler in setting up the "herniation" test cases (Rachmeler et al. 2007). The autoload PDL routine twisted_loop() generates a string declaration suitable for parsing with str2world, that describes a semi-elliptical loop with rectangular cross section and a prescribed total amount of twist. The loops are not constant-alpha; rather, they are quasi-Gold-Hoyle loops: the winding number about the central axis is the same for each fluxon in the loop.
Shot topologies may be found by shooting field lines through an analytic or other field solution from a boundary. Placing the initial field lines is key to getting the topology right.
The problem of placing fluxons at the boundary, given a measured or computed magnetogram, is equivalent to the mathematical problem of dithering. Several routines exist for placing the fluxons given a flux distribution at the boundary. A nice routine by Kankelborg (written in IDL) uses a divide-and-conquer strategy to place fluxons on a plane with an imposed magnetogram.
Two separate dithering codes exist in FLUX:
fluxon_placement_fs(): Floyd-Steinberg dithering places fluxons based on an error-diffusion algorithm that works along boustrophedonic rasters. fluxon_placement_hilbert(): Uses a linear error-diffusion algorithm that only propagates error along a single path through the image; but the path is a Hilbert curve, an area-filling fractal with a nice compromise between localization and no preferred direction.
Floyd-Steinberg is more widely trusted in the imaging community, and seems to work quite well for rough fields. But for smooth fields such as a priori models, the Hilbert placement does quite a bit better than the Floyd-Steinberg technique.
Two types of boundary are implemented: spatial boundaries that affect the whole arena and that can be can be free, line-tied (photospheric), or open; and fluxon end conditions that affect fluxon endpoints and can be line-tied, open, or plasmoid. The fluxon end conditions are part of model bookkeeping and are discussed at the end of this section.
Most solar simulations will use a planar photosphere or a spherical photosphere, with an outer boundary that is either free or open.
Below we describe the different kinds of boundary condition you can set up, and how to do it. You will need to understand how FLUX represents its objects -- for Perl scripts, they are represented as Perl hashes that contain specific keyword/value pairs, and this is the main way we describe them. For C programs, they are represented as C structures. Browse the Perl module sections of the manual for more information.
A free boundary is simply a lack of boundary. Unless you specify a boundary condition, the field value is extrapolated smoothly to zero at large distances. Fluxons are free to expand into the entire unoccupied space. To achieve a free boundary condition, do nothing.
Photospheric conditions are line-tied. Fluxons cannot cross through the photosphere. Fluxons may end at the photosphere, but only at line-tied flux concentrations that happen to be located on the photosphere. The no-crossing rule is enforced by the method of images and by explicit forcing. Normally, each fluxel (fluxon segment) interacts not only with real neighbor segments in the rest of the arena, but also with an image segment reflected through the photosphere (or a tangent plane of the photosphere). The image segment prevents the flux from the fluxel from expanding through the photospheric boundary, so that magnetic pressure prevents the fluxel from drifting down to the boundary itself. Because in extreme cases the discrete relaxation steps can take a vertex inadvertently (and disastrously!) through the photosphere, the code detects such steps and prevents the vertex from stepping all the way through the photoshere.
Photospheres can be planar, spherical, or cylindrical. The spherical and cylindrical cases are implemented by identifying the tangent plane closest to the subject fluxel, and reflecting through that tangent plane rather than through the analytic curve of the photosphere.
You can use up to two photospheric boundaries in a FLUX simulation. You set them up (in Perl) by defining the photosphere and photoshere2 fields of a Flux::World object, or (in C) by manipulating the corresponding PHOTOSPHERE structures contained in a WORLD structure. A PHOTOSPHERE contains a PLANE structure that can describe a plane in 3-space, and a type code. The PLANE contains six numbers that are used differently depending on the type of photosphere.
In the Perl side, a photosphere appears as a Perl hash with three keys:
origin: a 3-PDL containing the origin of the shape describing the photosphere normal: a 3-PDL containing the normal vector of the plane, or something else as described below type: a scalar containing a type code describing what the photosphere is.
There are currently four types of photospheric boundary:
- No boundary: set the type code to 0; all other values are ignored.
- Planar boundary: set the type code to 1; the origin and normal of the plane are specified as you might expect. Non-endpoint vertices are forced to the same side of the plane as the normal vector.
- Spherical boundary: set the type code to 2; the origin is the center of the sphere, and the 0 component of the normal vector is the radius of the sphere. The 1 component of the normal vector indicates whether the sphere is an inner boundary (the usual case, with the simulation outside the sphere; set the 1 component to be nonnegative) or an outer boundary (ie an enclosing ball, with the simulation inside; set 1 component is negative).
- Infinite cylindrical boundary: set the type code to 3; the origin is the center of the cylinder, and the normal vector points along the cylinder's axis. The length of the normal vector indicates the cylinder's radius. There are no end-caps on the cylinder. Only the interior of the cylinder can be used.
The boundaries are enforced at each relaxation step in two ways: compliant vertices interact with image vertices reflected through the boundary; and any non-compliant vertices are forced to be just inside the boundary at the end of the relaxation step. Vertices can become non-compliant through quirks of the relaxation, or through interaction with a different boundary -- for example, an open boundary (below) can move a vertex through a photosphere, making it non-compliant. At the end of the relaxation step the vertex is moved to the closest boundary-compliant part of space.
Open boundary conditions represent the transition to the solar wind, where outward flow is faster than the Alfven speed and field motions outside the boundary cannot affect the interior of the simulation. The open boundary is a large sphere. Fluxons that would pass outside of the boundary can be forced (by cutting and splicing if necessary) to end on the open boundary. The endpoint is allowed to evolve under the influence of the pressure force, but at each relaxation step it is forced onto the sphere by projection along a radial line.
Open boundaries are described using a special fluxon end condition, and two special flux concentrations are provided to receive open fluxons. On the Perl side, these flux concentrations are the fc_ob and fc_oe fields of the Flux::World object, and you can open a fluxon simply by attaching it to one of these flux concentrations. They have negative ID numbers to distinguish them from normal flux concentrations. Rather than explicitly reattaching fluxons to the open flux concentrations mid-simulation, you can use an "auto-open" feature to open up fluxons that pass through the open fluxon boundary.
To set up open boundary conditions for your simulation, set the auto_open flag in the Flux::World, and make sure that the special open flux concentrations (accessible as the fc_ob and fc_oe fields of your Flux::World arena) have the right location and radius. (The open-flux surface is defined by a sphere centered on the location of the open flux concentrations, with radius set by the locale_radius of those concentrations. The two special open concentrations should have the same locale_radius and x, so if you change one you should change the other.) Fluxons that venture outside of the open sphere will automatically be cut and opened, and can form U-loops that will propagate outward through the boundary.
You can also do special tricks by defining your own open flux concentrations -- an open flux concentration is any concentration with the "fl_b_open" end condition in its bound field, so you can create one explicitly if you choose.
To eliminate the open boundary condition, set the auto_open flag to 0 and make sure none of your fluxons terminate on an open flux concentration.
Plasmoids are "loops" of magnetic force that are not explicitly connected to any boundary (though they may be held in place by geometric linking to other fluxons). They are implemented using the "plasmoid" fluxon end condition (see below). Plasmoids can be created explicitly or by reconnection of a fluxon to itself. Two special flux concentrations (accessible as the fc_pb and fc_pe fields of a Flux::World) are intended to hold all the plasmoid fluxons in the entire simulation, though you can define your own plasmoid concentrations for special tricks or accounting.
To make a plasmoid, attach the beginning of a fluxon to the fc_pb concentration and the end of it to the fc_pe concentration. You should also set the fluxon's plasmoid field to unity.
You can tweak the way that endpoints of fluxons are handled, by changing the bound field of any Flux::Concentration object. The default end condition for new flux concentrations is set in the default_bound field of the Flux::World, and existing concentrations can be modified explicitly through their bound field: $world->concentration($n)->{bound}="boundary_name" or just $conc->{bound}="boundary_name". (You can get a list of the valid boundary condition names by submitting an erroneous one: $conc->{bound} = "foo" will generate a suitable error message, assuming that $conc contains a Flux::Concentration.
In Perl, the end conditions are represented as specific strings that name the C subroutine used to implement the end condition. They are:
-
fl_b_tied_inject: ties the endpoint of the fluxon to the location of the flux concentration; injects new vertices if the fluxon stretches so that the first non-endpoint vertex is outside the locale_radius of the flux concentration.
-
fl_b_tied_force: ties the endpoint of the fluxon to the location of the flux concentration; forces the first non-endpoint vertex to lie on a sphere centered on the flux concentration, with the locale_radius of the concentration.
-
fl_b_open: causes a flux concentration to be associated with an open boundary condition -- the endpoint is forced to the surface of a sphere with the locale_radius of the concentration, centered on the concentration itself.
-
fl_b_plasmoid: causes a flux concentration to hold a plasmoid with no explicit anchor to the rest of the simulation (except for geometric linking). The endpoint locations are copied from the second-to-last vertex locations on the other end of the fluxon, forcing the fluxon to be a free-floating plasmoid. If one endpoint of a fluxon is connected to a plasmoid concentration, the other end should also be connected to a different (opposite sign) plasmoid concentration. The World comes with two plasmoid concentrations explicitly defined, in the
fc_pbandfc_pefields.
FLUX supports a collection of modular physical forces, to modify the physical system being simulated. You must select a set of force laws at run time. From the Perl front end, you can do that by calling the forces method, as in $world->forces( \@forcelist ), or (preferably) using the hash interface: $world->{forces}=\@forcelist. \@forcelist is a list of names of force laws. Specific named force laws are declared in physics.c and its header file physics.h. You can print a current list from the code itself by assigning an erroneous one: $world->{forces} = ["foo"] will throw an error message that includes all the current valid names.
To use FLUX to find nonlinear force-free fields, you need a minimum of three forces: a magnetic pressure force, a magnetic curvature force (those first two are just components of the Lorentz force), and a vertex pseudoforce. For force-free field extrapolation, you can use the original field-normalized "forces", which calculate force divided by magnetic field strength.
The current recommended set of force laws for force-free field calculation are:
- f_pressure_equi2b - field-normalized magnetic pressure gradient force
- f_curvature - field-normalized curvature force
- f_vertex4 - field-normalized vertex distribution pseudoforce
Several force laws are available -- some of them (such as "f_pressure_equi") are previous versions of currently recommended laws, others (such as "b_eqa") actually calculate other quantities than force, and still others (e.g. "f_p_eqa") are true forces rather than field-normalized forces. The field-normalized and true forces should not be mixed in a single relaxation.
Relaxation requires iteratively calculating net forces and taking small steps in the direction of those forces. Relaxation is managed via simple_relaxer, a Perl/PDL subroutine that makes calls to the C relaxation library.
The current setup seems to work well with a dtau setting of 0.25. This yields relatively rapid convergence and also good stability. dtau can be set as high as 0.5 in some circumstances.
You can maintain close to the ideal number of vertices by periodically calling fix_curvature (which is accessible from within Perl). Fix_curvature takes a maximum curvature and a minimum curvature to allow, and inserts and removes vertices to keep the curvature at each vertex within those values.
Calling fix_curvature every 100 or so steps during the initial stages of relaxation, with values of (0.3, 0.1), seems to work quite well. That limits the per-vertex curvature to be between about 6 and about 20 degrees, and speeds convergence considerably compared to just adding a large fixed number of vertices at the beginning.
A word of warning -- instabilities that might be tolerable without periodic curvature fixation tend to grow out of control with frequent curvature fixing! If you're having problems with convergence, try backing off on the fix_curvature, or leaving it out altogether.
FLUX is a relaxation code, so the fundamental operation is a series of relaxation steps that successively approach an equilibrium by moving the vertices of the simulation. The forces driving the vertices are tuned so that when the forces are at equilbrium, the system is a good approximation of the physical state.
At each step, all the individual force terms acting on each ith vertex are summed to produce a total force:
The vertices are relaxed in the direction of the force:
where α is a proportionality factor and Δτ is a relaxation step size.
Simple force-proportional stepping doesn't work efficiently or even adequately for a system like FLUX: field lines are far apart where forces are weak, and close where forces are strong -- so that the relaxation step Δτ would have to be tiny to prevent instability in strong-field regions, and simultaneously huge to prevent taking forever in weak-field regions. FLUX overcomes that by multiplying the step size by locally calculated parameters. Each of several coefficients is raised to a power, and they are all multiplied together to create α:
where the ξk are several locally calculable quantities and the γk are adjustable power laws intended to make the actual step size scale like the needed step size.
The power laws need to be selected carefully to avoid instability while simultaneously not sinking the simulation into molasses. They are accessible via the Perl hash interface. The current scaling coefficient powers are:
- d (hash element:
scale_d_power) - the distance from the current vertex to the closest other vertex. Useful for accelerating weak-field motions and slowing strong-field motions...
- b (hash element:
scale_b_power) - the local magnetic field strength. Related to d -- b should be something like d^-2 -- but useful.
- s (hash element:
scale_s_power) - the stiffness coefficient: the ratio of the magnitude of the resultant force on the vertex, to the sum-of-magnitudes of all the component forces on that vertex. S is a measure of how relaxed the current system is. Useful for preventing instabilities near final relaxation.
- ds (hash element:
scale_ds_power) - the similarity of the force on the current vertex to those of the adjacent vertices on the same fluxon. Calculated as the ratio of the magnitudes of the sum and difference of the local and nearby forces. Useful for accelerating relaxation of fluxons with many vertices, where each vertex may have nearly the same force on it as its neighbors.
Currently recommended values for field-normalized relaxations: scale_s_power=0, scale_d_power=2, scale_b_power=0, scale_ds_power=0 (but not more than 0.25).
All relaxation codes have a specific efficency problem: near equilibrium, the force law can be described as a linear expansion around the equilibrium: with n free parameters, the force law Jacobian is an n×n matrix. The eigenvectors of the matrix that have high eigenvalues converge rapidly and require small relaxation steps to avoid instability. The eigenvectors that have low eigenvalues converge slowly. The ratio of largest to smallest eigenvalue sets the number of relaxation steps required to converge. In general, those eigenvalues can range over many orders of magnitude.
Because FLUX uses local gradients to generate forces, the lowest eigenvalues correspond to broad spatial displacements: motions of many vertices in the same direction. FLUX includes a mechanism to boost those eigenvalues in the relaxation. At each relaxation step, the "naive step" for each vertex is calculated using the force laws as described, and stored in the vertex structure. Then a second pass is made through the data, to accumulate local average motions. The final step made by each fluxon is the sum of several spatial averages of calculated naive steps over successively broader areas.
The averaging takes place by "neighbor generations" - the 0th generation consists of the vertex itself, and the 1st generation consists of the vertex's neighbors. The 2nd generation is neighbors of neighbors (that are not themselves the original vertex or one of its neighbors). Each generation's average displacement is calculated separately, and the final step is the sum of these averages, each scaled by an independent coefficient. This raises the rate at which large spatial scale eigenvectors relax.
The coefficients are set by the Perl/PDL "coeffs" hash field in the Flux::World object's tied hash. You assign an array of coefficients, and the code iterates out to the size of the array that you specify.
Currently, neighborhood acceleration is experimental. To switch it off, set the "coeffs" field to [1].
Reconnection is handled via a thresholding technique. Several "reconnection criterion" routines are available, that check for certain local physical heuristics to exceed a threshold between each pair of vertices. If the threshold is exceeded, the vertex pair is marked for reconnection and, er, reconnected.
Reconnection handling is generally independent of relaxation: it is recommended to relax, then check for reconnection (using the routine global_recon_check)
(main page: User Manual: Code)
FLUX is a scaled Eulerian (step-style) relaxer. At each relaxation step, the code calculates the net force on each vertex, and then moves each vertex in the direction of the total force on it. The forces are calculated according to a set of parameterized force laws called out in the run specification. The geometry relaxes (hopefully) toward a force-free state.
There are two parts to the code: a library written in C and a control interface into Perl/PDL. Some IDL interface code has also been written.
FLUX uses dynamically allocated structures rather than regular arrays to handle its data. The structures are quasi-objects, although the code is written in C (and not C++ nor Objective C). The C library handles the most computing-intensive parts of the simulation, using its dynamic data structures. The PDL interface code allows access to the data structures as tied perl hashes: although the underlying data is contained in a linked collection of C structures, they are accessible as a linked collection of Perl hash refs.
Main article: User_Manual:_Data_Structures.
Data within FLUX are represented as a collection of dynamically allocated structures. There is a "glue layer" of code that ties these structures to Perl hashes, so you can manipulate (most aspects of) your FLUX simulation directly from a Perl script. The main object types are WORLD (which represents the simulation as a whole), FLUX_CONCENTRATION (which represents a fluxon endpoint), FLUXON (which represents a fluxon), and VERTEX (which represents a single fluxel).
Several other data structures represent geometric objects such as 3-D and 2-D vectors, hull points, etc.
There are no permanent global variables in the C code. (Some global variables exist as placeholders for callback routines).
There are five basic parts to the C library:
data.c: contains low-level interface routines for interacting with the data structures
io.c: contains a simple set of routines for generating and parsing ASCII files that describe FLUX arenas.
geometry.c: contains a low-level 3-D geometry library including quite complex tasks such as the non-Cartesian projection and hull algorithms.
physics.c: contains all of the physics-specific code: forces, pseudo-forces, and magnetic field calculation.
model.c: contains most of the entry points necessary to run a FLUX model, notably the global_relax_step() function that calculates forces and updates positions of all vertices.
(Main article: User Manual: Code: model.c)
The model infrastructure segment of the library handles tree walking, step calculation, neighbor candidate grouping, and the other bookkeeping aspects of the simulation.
Relaxation framework
A typical relaxation step works as follows:
- Loop over all vertices in order, and for each one:
- Expand the current neighbor list (by local dilation or, if selected for "global" checking, to the whole simulation).
- Add any necessary pseudovertices (from photospheric boundaries, etc.) to the current vertex's neighbor candidate
- Find the current neighbors by projecting into the vertex's following segment's perpendicular plane, and carrying out the hull algorithm in physics.c.
- Call all of the force routines, in order, on the prepared vertex
- Use the force routines' output to estimate the force on the vertex
- Muliply the force by the standard coefficients raised to the global scaling powers, to get a planned step.
- Make a second pass through all the vertices, and for each one:
- Loop over successive ranks of neighbors and:
- Find the mean planned step for the current rank, multiply by the expansion coefficient for the current rank, and add to a step accumulator
- Shift the vertex by the global delta-tau times the accumulated step.
- If the vertex at the end of a fluxon, apply fluxon-end boundary conditions as required.
- Loop over successive ranks of neighbors and:
The first pass finds the relaxation force and estimates how the vertices should move toward relaxation; the second pass is used to accelerate relaxation of large systems by increasing the effective grid speed. For non-accelerated relaxations, there is only one coefficient (for the zeroth tier of neighbors), and the step is proportional to the locally calculated forces.
Additional operations such as rendering, reconnection, etc. are generally applied individually from the Perl side.
Boundary conditions
Boundary conditions come in two types: global boundaries (such as a photosphere or globally defined open surface) and fluxon-specific boundaries (such as line-tying, open-surface tying, or plasmoid conditions).
The global boundaries such as the photosphere are implemented by the method of images: each fluxel gets, as a candidate neighbor, an image of itself reflected through the boundary. If the image is close enough to be a true neighbor, the fluxel interacts with the boundary; otherwise, it doesn't.
Fluxon-specific boundaries are implemented as a final correction to the position of the end vertices in each fluxon: line-tied conditions force the first and last vertices to specific positions; open conditions force the last vertex onto the surface of a large sphere; and plasmoids force the first and last segments to overlap by copying the second-to-last vertex's position and the second vertex's position into the first and last vertices, respectively. Global specification of the open surface is achieved by forcing all open fluxons to begin or end (as appropriate) at specific flux concentrations that are part of the global WORLD data structure.
For more detail, see User Manual: Code: model.c.
(main article: User Manual: Code: geometry.c)
The geometry package provides two main functions: it contains the hull algorithm used to find the Voronoi cell of each vertex for the magnetic field and force-balance calculation; and it contains an extensive set of utilities for manipulating 2- and 3-vectors. That is important because it aids writing the physics routines in a co-ordinate free manner. While vectors are stored in Cartesian coordinates, the higher level portions of the code (such as the hull calculation and the physics modules) are written in a coordinate free manner using vector primitives implemented in geometry.c.
The package is minimal but implements norm, scalar multiplication, addition, subtraction, dot product, cross product, matrix multiplication and construction of a particular family of projection matrices. It is intended to encapsulate all coordinate-dependent operations, so if you find yourself using the Cartesian components of a vector (say, in a new physics.c module) then you should isolate those operations and put them into geometry.c.
For a full list of entry points and functionality, you should read the main article on geometry.c.
(main article: User Manual: Code: physics.c)
The physics package contains the microphysics of the simulation. It is a collection of "physics modules" that can be called to act on the force accumulators and other data in each vertex, in the innermost loop of the simulation. Physics modules have a standard interface (defined in physics.h) and are called in order from a queue that is set at runtime. From PDL, the Flux::World::forces method sets and/or retrieves the current force list of a world, as human-readable ASCII strings. (The internal storage format is as a set of function pointers).
For a description of what forces are currently available, and the recommended set, see User Manual: Code: physics.c.
The main interface for using the C library is a set of Perl objects that rely on the PDL extension to Perl. Perl incorporates object oriented extensions, so that it is possible to define new "objects" of specific types, and "methods" (type-specific subroutines) to manipulate them.
The main data structures of the C library -- WORLD, FLUX_CONCENTRATION, FLUXON, and VERTEX -- are implemented as specific Perl object types -- Flux::World, Flux::Concentration, Flux::Fluxon, and Flux::Vertex, respectively -- that can be manipulated from within Perl scripts. Flux objects are themselves Perl tied hashes, so that Perl scripts can access the elements of each data structure as elements of a hash as well as using the supplied access methods. Thus, Perl control scripts have (nearly) full access to the tree and neighbor structures in a Flux::World.
In addition to the object definitions, Flux provides a set of PDL autoload scripts (analogous to .pro files in IDL) that provide functionality at the script level.
Individual articles discuss in detail the functionality of each object type, and of the autoload library. They are:
Flux.pm: a few global utilities of use to all the Flux objects, including the hash-tying mechanism
Core.pm: A dymamic linker that services the other modules
World.pm: the Flux::World object, corresponding to the WORLD data structure in the C library
Fluxon.pm: The Flux::Fluxon object, corresponding to the FLUXON data structure in the C library
Vertex.pm: The Flux::Vertex object, corresponding to the VERTEX data structure in the C library
Concentration.pm: The Flux::Concentration object, corresponding to the FLUX_CONCENTRATION data structure in the C library
Autoload files: A collection of useful subroutines, stored in a PDL autoload directory for simple access
A discussion of the energy calculation can be found on the Energy Calculation page.
More on running the code soon...
Worked Examples may be found here.