User Manual: Code: physics.c

From FLUX

Jump to: navigation, search

(Return to Code or the User Manual)

physics.c contains the "meat" of the simulation, in the form of standardized, modular local physics calculations. Each of the routines in physics.c may be linked into a world's "force list", and all of the force list is executed in order on each VERTEX after the hull code is executed. Each physics routine performs some relevant calculation, and adds forces to the various accumulators in the VERTEX. The routines have full access to the VERTEX, and can do other things than force calculation -- for example, some routines estimate the B field and store it.

Recommended physics routines

The currently recommended combination for force-free field calculation is:

f_pressure_equi2b
f_curvature
f_vertex5

Physics code calling environment

Each physics routine is described in a static array, FLUX_FORCES, at the top of physics.c. FLUX_FORCES lists the name of the routine, together with a brief description and a function pointer. It's used to translate force name (in the Perl layer) to function pointer (in the forces list in the WORLD structure), and vice versa.

The individual routines should be declared like so:

void x_name(VERTEX *v, HULL_VERTEX *verts)

where the x denotes what the thing is calculating (so far, 'f' for force, 'e' for energy, or 'b' for magnetic field), and the name is a mnemonic.

When called, the verts array contains the hull vertices for the neighbors listed in the VERTEX (and associated with the line segment following the vertex).

The force routines load the force per unit length of segment (not true force) into the appropriate force accumulators

Vertices and geometry

Forces should work in the perpendicular plane. On call, the the v->a and v->r fields contain the 2-D radial coordinates of each neighbor, and the v->scr field contains its 2-D Cartesian coordinates.

The verts hull vertex structures contain information about the boundary associated with the corresponding neighbor, and the actual vertex (ordinary or open) following the associated neighbor. The a_l and a_r elements indicate the absolute 2-D angle of the vertex when seen from the right (as in the current neighbor) or from the left (as in the next neighbor). Those values are the same for ordinary vertices, and different for open vertices (for which there is no intersection because the adjacent Voronoi cell walls diverge rather than converging).

Forces

Each VERTEX has several force accumulators of interest to force routines. After calculating your force, you should add it to the appropriate force and magnitude accumulators. Remember that the quantity involved is force per unit length, not force itself.

f_s: vector accumulator for forces that occur at segment centers (halfway between vertices), such as pressure forces and gravity.
f_v: vector accumulator for forces that occur at vertices, such as the curvature force.
f_s_tot: scalar accumulator for segment force magnitudes
f_v_tot: scalar accumulator for vertex force magnitudes

In addition each vertex has fields to hold magnetic field calculations:

b_vec: the calculated vector magnetic field

b_mag: the magnitude of the magnetic field

energy: the magnetic field energy associated with this segment

Future force laws may involve plasma forces, which include parameters like temperature and density of material, and velocity along the field line. Those quantities wil be stored in the vertex but are not yet (in version 1).

Lists of force routines

Because FLUX is development code, many force routines have been written for it, not all of which are useful. In particular, at least two sets of magnetic force routines exist. The first set, with long names like f_pressure_equi, return force per unit length (not force) divided by the local magnetic field strength. Those force laws are the best validated and are used for the force free field calculations. They are not suitable for generalization to systems where plasma forces are important. The second set, with short names like f_p_eqa, are genuine forces per unit length (not just force). The two types of force law should not be mixed in a single simulation!

B-normalized force routines

The currently (25-Sep-2007) recommended set of force laws is f_curvature, f_pressure_equi2b, and f_vertex4. These laws need to have the d scaling power set to 2, all others to zero.

f_curvature

(Recommended)

Normalized curvature force, located at the vertex. The force per unit length per unit magnetic field is $\left(\hat B\cdot\bold{\nabla}\right)\hat B$ ; that is a curvature (inverse radius). That's just the bend at the vertex, divided by the associated segment length.

f_curvature2

Modified normalized curvature force has an extra fudge factor to avoid runaway when the small-angle approximation breaks. Tends to cause instability instead, in extreme cases.

f_pressure_equi

This is $-\nabla_\perp\left(|B|\right)/|B|$ . It uses the approximation:

$\nabla_\perp\left(|B|\right)/|B| \sim \sum_i\left(\hat n_i \Delta\phi_i\right) / \left(\pi r_i\right)$ where $\hat n_i$ is the (2-D) unit vector normal to the Voronoi cell boundary segment, $Δφ$ is the range of angles subtended by the Voronoi cell boundary segment, and $r i$ is the radius to the segment.

This approximation (the Kankelborg/DeForest discretization) underestimates the magnetic pressure gradient by approximately a factor of two and is not recommended. See below for more.

f_pressure_equi2

The formulation in f_pressure_equi estimates one of two terms to the gradient of B in terms of the density of field lines -- the variation in fluxon radial density with respect to radius. The other, variation in fluxon angular density (suitably normalized) with respect to radius, is ignored. In equilibrium, most cells are close to circular and hence the two terms are approximately equal. So we estimate the full gradient by simply doubling the first term.

This force law overestimates the pressure force on the last field line, which has an open Voronoi cell. That causes magnetic systems to grow much larger than they should.

f_pressure_equi2a

This formulation is similar to f_pressure_equi but uses a slightly different approximation. f_pressure_equi treats the magnetic pressure as acting normal to the Voronoi cell walls at every location. f_pressure_equi2 treats it as acting radially with respect to the central fluxon, by including force components parallel to the Voronoi cell wall. For closed Voronoi cells this is a difference that makes no difference: the parallel forces can be shown to cancel around the cell. For open Voronoi cells, this reduces the outward pressure gradient, effectively allowing the outside field lines to press inward against the central fluxon.

This force law seems to approximate the field forces on the last field line, but underestimates the pressure force on the last field line by a small amount -- enough to shrink our test Low & Lou solution by about 20% upon relaxation.

f_pressure_equi2b

(Recommended)

This formulation is similar to f_pressure_equi2a, but does not double the parallel forces as it does the perpendicular forces on the Voronoi cell walls. That doesn't affect closed cells but slightly increases the outward pressure gradient on the last field line in free-boundary situations.

f_vertex

Vertex distribution force includes a vertex repulsive force, a curvature attactive force, and a proximity-to-other-fluxons attractive force, in an attempt to evenly distribute vertices. Appears to work reasonably well but fails for two reasons: (1) it is grossly mismatched from the vertex placement algorithm, so that long time relaxation with occasional curvature fixing yields a net vertex flow; and (2) the force law doesn't scale properly with distance, so that very large flux systems become unstable.

f_vertex2

Do not use; kept for historical reasons.

f_vertex3

Do not use; kept for historical reasons.

f_vertex4

This is a two-component vertex distribution force that appears to scale the same way as the magnetic forces for the B-normalized case. It includes a 1/r^2 repulsive force between adjacent vertices (normalized to the square of the vertex length, so that the resultant force per unit length is independent of physical size). The curvature attractive force per unit length is also length independent, being determined only by the angles subtended by the vertices at either side of the subject vertex.

f_vertex4 tends to distribute the vertices too regularly, as it neglects the neighbor distance entirely. That leads to stability problems near concentrated footpoints in the photosphere, which propagate out to form the "seaweed effect" as the overlying field structures sway under the influence of instabilities at the footpoints.

f_vertex5

(Recommended)

This is a two-component vertex distribution force that works mainly by proximity to other fluxons: the main component is dependent on the ratio of segment length to nearest-neighbor distance, so that vertices tend to distribute themselves to equalize this ratio. There is an additional perturbation that is proportional to 1/r^2 along the fluxon (but scaled by 0.1), to help regularize the vertex behavior. f_vertex5 is the currently recommended vertex distribution force.

Non-B-normalized force routines

Non-force physics routines

Reconnection

Reconnection is handled using a similar mechanism to force calculation. Reconnection routines accept a VERTEX pointer and a collection of parameters that are normally stored in the WORLD object. They are called by the routine vertex_recon_check, which in turn is called by the main entry point global_recon_check/code>. Both of those are in model.c.

Currently there is only one reconnection criterion; more are planned.


[edit]
rc_a_ad2
This is a reconnection criterion that causes reconnection between adjacent segments if two criteria are met: a minimum angle and a minum angle/distance^2.  Both angles are specified in radians.

[edit]
rc_a_ad2_h
Same as rc_a_ad2, but with a sliding current threshold that falls off exponentially with altitude; this simulates a static atmosphere in a plane-parallel case

[edit]
rc_a_ad2_loc

Same as rc_a_ad2, but with a sliding current threshold that falls off near a particular location, to induce reconnection only in a particular locus.

Retrieved from "http://flux.boulder.swri.edu/wiki/index.php/User_Manual:_Code:_physics.c"

Views

Personal tools

Search

Toolbox