Frequently Asked Questions

Is it possible to display fluid markers in VisIt?

Yes, with the minor caveat that these markers need to be evolved inside a simulation. Fluid markers are managed by the IBTK::MarkerPatchHierarchy class and can store data in H5Part files (a common marker point format which VisIt knows how to read). IBFE example 4 demonstrates how markers work.

What do the different input files do?

For an IB model, in addition to the input2d/input3d files, there are a number of files that specify the initial conditions and material properties of the immersed elastic structures. The only such file that is required is the “.vertex” file, which specifies the initial positions of the IB points.

There are several additional files that may be optionally specified:

A “.spring” file can be used to specify an essentially arbitrary network of linear or nonlinear springs that connect the various IB points. Each “spring” connects precisely two IB points.

A “.beam” file can be used to specify a collection of linear beams with a preferred curvature. Each “beam” connects precisely three IB points.

A “.target” file is used to specify IB points that are, by default, “tethered” by stiff linear springs to their initial positions. Note that it is also possible to update the location to which each target point is tethered to within a simulation, so as to provide an approximately-prescribed motion for the target points. (See below.)

A “.anchor” file is used to specify “anchor” points. In IBAMR, anchor points are IB points that neither apply force to the fluid nor interpolate velocity from the fluid. They are literally anchored in place. (Note: Anchor points are distinct from target points. Target points can be thought of as being tethered to anchor points.)

A “.mass” file is used to specify any additional mass associated with the IB points. For such files to have any effect, it is necessary that the IB solver be run in “penalty-IB” mode.

How can I have time-dependent target point locations?

Take a look at https://ibamr.googlecode.com/svn/examples/moving_plates. This example is for the specific case of time-dependent target point locations. The linked code has a routine (update_target_point_positions) that loops over the local IB points (i.e., the IB points that are “local to the processor”). The pseudo-code is roughly:

for each local IB point
{
    check to see if the point has a target point spec
    if yes: update the material properties w/that target spec
}

The properties for the target points are stored in IBTargetPointForceSpec objects. Functionality has been added to the most recent IBAMR versions to simplify the code required to access these objects. In particular, for each Lagrangian (IB) point, you can do:

Pointer<T> spec = node_idx.getStashData<T>();

to get a pointer to any object of type T that is associated with the Lagrangian point.

How can I have time-dependent spring constants or resting lengths?

If you use class IBStandardInitializer to initialize the Lagrangian structure and provide a “.spring” file, then the first point in each spring will have associated with it an object of type IBSpringForceSpec. You can access that object via:

Pointer<IBSpringForceSpec> spring_spec = node_idx.getStashData<IBSpringForceSpec>();

Note: Each spring is associated with ONLY ONE of the two ends of the spring. By default, this is the “master” point index, which is simply the lower of the two indices. This means that it is possible for some IB points not to have any force spec objects, even though those points might be connected to another point by a spring (or a beam). If there are not any objects of a specified type T associated with a particular node, then the function getStashData<T>() will return a null pointer. You should always check to see if the returned pointer is non-null before trying to do anything with it.

To update the properties of the objects that describe the springs will require a function similar to the update_target_point_positions function mentioned in the previous Q/A. Refer to the documentation for class IBSpringForceSpec to see what properties you can access and modify.

Which part of the IBAMR takes the most runtime?

Cartesian grid linear solves generally dominate the overall run time, but this can be highly application dependent. Improvement in linear solver performance would have a big impact on overall code performance.

What method is used for adaptive mesh refinement?

The basic methodology is called “block-structured adaptive mesh refinement.” Block-structured AMR is essentially a description of the data structure used to describe the locally-refined discretization. Block-structured AMR grids are comprised of hierarchies of nested levels, each of which are comprised of one or more rectangular grid patches. Most of the algorithms and data structures related to AMR are provided by the SAMRAI library (https://computation.llnl.gov/casc/SAMRAI).

AMR grids are generated by SAMRAI using the Berger-Rigoutsos box generation algorithm.

What is the numerical method to solve the incompressible Navier-Stokes equations?

There are currently two different incompressible flow solvers implemented in IBAMR.

One is a cell-centered approximate projection method-based solver that uses an explicit second-order Godunov scheme to handle the nonlinear advection terms that appear in the momentum equation. The current version of IBAMR always uses the xsPPM7 variant of the piecewise parabolic method (PPM) for the Godunov scheme. The basic timestepping methodology is a predictor-corrector scheme that is a combination of the explicit midpoint rule (for the advection terms), the explicit trapezoidal rule (for the IB forcing terms and IB advection equation), and the implicit trapezoidal rule (for the viscous terms and the incompressibility constraint).

The other solver is a staggered-grid scheme that uses a projection method-based preconditioner for the incompressible Stokes equations. Again, xsPPM7 is used to approximate the nonlinear advection terms. The basic timestepping methodology is truncated fixed-point (Picard) iteration that uses the explicit midpoint rule for the advection, IB forcing, and IB advection equations, and the implicit midpoint rule for the viscous terms.

Both solvers are formally second-order accurate (i.e., second-order accurate in the discrete L1 or L2 norms), but suffer from a local reductions in accuracy at coarse-fine interfaces in the AMR grid (i.e., are only first-order accurate in the maximum norm). In our experience, such reductions in accuracy are readily-apparent only at low Reynolds numbers. At higher Reynolds numbers, both schemes are effectively second-order accurate.

The staggered-grid scheme is recommended for two reasons: (1) Its handling of physical boundary conditions is more accurate and robust, and (2) it possesses dramatically improved grid-independence and volume-conservation properties when used with the IB method.

Are there different solver options that should be used for low Reynolds number flows as compared to high Reynolds number flows?

To obtain good performance at low Reynolds numbers, it is necessary to use a FAC/multigrid preconditioner for the implicit discretization of the viscous terms.

What grid refinement ratios should be used with IBAMR?

In principle, for the cell-centered solver, you can set the refinement ratio between AMR grid levels to be any even number. (Odd refinement ratios used to work correctly but have not been tested in several years.) Note that you can also use different refinement ratios between different levels of the grid.

However, the discretizations that are implemented within IBAMR have a better formal order of accuracy for even refinement ratios of at least 4.

Parts of the staggered-grid AMR discretization implemented in IBAMR require that the refinement ratio be a power of 2. As in the cell-centered case, refinement ratios of 4 or larger yield a higher formal order of accuracy than a refinement ratio of 2. In particular, with a refinement ratio of 2, there are O(1) errors at coarse-fine interfaces in the viscous terms and O(h) errors at coarse-fine interfaces in the pressure gradient. With a refinement ratio of 4, there are O(h) errors at coarse-fine interfaces in the viscous terms and O(h^2) errors in the pressure gradient.

How do periodic boundary conditions work?

Periodic boundary conditions are handled by SAMRAI. If periodic boundary conditions are enabled in any direction (via the “periodic_dimension” flag), they take precedence over all other boundary condition-related settings. More precisely, when periodic boundary conditions are enabled, all “physical boundary condition”-related settings are ignored.

SAMRAI, which provides many of the basic AMR algorithms and data structures used by IBAMR, divides patch boundaries into three types: (1) physical boundaries, (2) periodic boundaries, or (3) internal boundaries. The physical boundary condition-related code is only called for patches that touch the “physical” boundary.

To reiterate: If periodic boundaries are enabled, then the physical boundary condition routines are never called, and any prescribed boundary values are ignored.

Do I compute the Lagrangian force or the Lagrangian force density?

Although force densities are more natural in the continuum setting, the Lagrangian data structures implemented in IBAMR are not sophisticated enough to use force densities. To simplify the software interface, we decided to make all Lagrangian forces actual forces (not force densities). Note that this convention affects the units required for spring constants, bending stiffnesses, etc.

Essentially, if you have a Lagrangian elasticity model that is phrased in terms of force densities, you must convert it into a model that is phrased in terms of forces in order to use this model with IBAMR. If F is the force passed to IBAMR and F~ is the original force density, what you need to compute is something like F(s,t) = F~(s,t)*dV, where dV is the appropriate (discrete) volume element.

We are in the process of adding support for finite element-based elasticity models that will allow elasticity models that are actually stated in terms of force densities.

Is it possible to include additional advected-and-diffused scalar (or vector) quantities in an IB model?

See the “navier_stokes” example located at http://ibamr.googlecode.com/svn/examples, which shows how to add an “auxiliary” advected quantity to the staggered-grid solver.

It is also possible to handle reaction-advection-diffusion equations. The interface to the advection-diffusion solver allows you to provide an optional source term that can be used to handle the case where it is sufficient to treat the reaction term explicitly.

Is it possible to obtain a non-symmetric solution for problems where one might expect to retain symmetry?

Unfortunately, yes. Probably the first place where loss-of-symmetry can appear is in the AMR grid generation algorithm. The grid generation algorithm implemented in SAMRAI does not guarantee that “symmetric” patches are generated for “symmetric” cell-tagging patterns. Non-symmetric grids are especially likely for cases in which there are “holes” in the tagging pattern.

Even if the grid retains symmetry, there are subtle ways in which symmetry can be lost. In particular, there are several implicit orderings in the Cartesian grid linear solvers that can result in a loss of symmetry in the numerical scheme. There is also the potential for loss of symmetry in apparently simple operations such as computing forces on the Lagrangian mesh and spreading forces from the Lagrangian mesh to the Cartesian grid!

We’ve done a number of symmetric test cases to try to catch any obvious symmetry-related bugs, and we suspect that “extreme” loss of symmetry is indicative of numerical under-resolution.

How should the timestep size be chosen?

The incompressible flow solvers are only stable up to a CFL number of 1. However, in general, the IB method imposes some additional timestep restrictions that are related to the explicit treatment of the Lagrangian forces. When the timestep size becomes too large, this generally results in a rapid onset of instability that culminates in a reduction of the timestep size in order to accommodate the prescribed CFL number.

Although these timestep restrictions result from the stiff springs and beams associated with the Lagrangian elasticity model rather than the CFL condition of the advection scheme, the easiest way to diagnose the onset of an instability is by monitoring for “CFL violations”, which will typically occur soon after the IB method goes unstable.

To “optimize” the timestep size, you can try to systematically double the timestep size in the input file, doing so until you get a CFL violation. By default, IBAMR solvers will terminate when a CFL violation occurs.

Why are there pressure initial conditions?

Pressure initial conditions are only for purposes of VisIt output — pressure initial conditions are not used in any way by any of the fluid solver.