Module multimaterial 

Coupled materials: several fluids on one grid

The crate merlin-multimaterial advances several materials on one grid together. Each material has a field over the whole grid, its own equation, scheme and boundary conditions, and a level-set function that says where it is. It is the part of phase 9 that moves the interface and keeps a state for each material on both sides of it. What happens at the interface (the matching of pressure and velocity, the ghost states built from the interface response of materials) is phase 9.3, whose first step is described below: here the ghost cells of a material are the extension of its own real state, which is the baseline 9.3 improves on. There is no problem-file syntax for several materials yet, so every case below is built in code (the tests, and the example crates/merlin-multimaterial/examples/bubble_in_flow.rs).

The crate depends on core, systems, finite-volume and the level set, and names no writer (crates). It is built to take solids later without a change of the solver: a material is a MaterialSlot, a trait object that names no fluid model and no pressure; ConservativeSlot takes any equation (one with a velocity, MaterialEquation) and any flux scheme.

One step

A MultiMaterialSolver owns the fields (the state of each material, then its level-set function), the clock and one Runge-Kutta loop over all of them (Shu and Osher 1988): every stage brings all the fields to the stage before any is evaluated at the next, so what one field reads from another is the value of the stage and not of the start of the step. With one material it is the run of Solver to the last bit. With several, a stage does the following.

1. Who owns what. The owner of a cell is the material whose function is largest there, a tie going to the lower index. The layer of a cell for a material is 0 where it owns the cell and otherwise the number of face-to-face steps to the nearest cell the material owns, up to the width of the band (the reach of the widest scheme plus two cells by default); beyond that it is one more than the width. The layers come from the owner map alone, one pass per layer with a halo exchange between the passes, so they do not depend on the blocks or the threads (LevelSets).

2. Ghost states. The band of each material, the cells with a layer from 1 to the width, is given a state from the material’s own real cells. ExtrapolateAlongNormal solves the extrapolation equation (Fedkiw, Aslam, Merriman and Osher 1999; Aslam 2004)

   q_\tau - \mathbf n\cdot\nabla q = 0,\qquad \mathbf n=\frac{\nabla\phi}{\lvert\nabla\phi\rvert},

   with the first-order upwind rule at a Courant number of one, which is its steady state, one layer at a time: a cell of layer $k$ takes the average of its neighbours of layer $k-1$ on the side of increasing $\phi$, weighted by the components of $\nabla\phi$. A cell reads only cells of lower layers, so the result depends neither on what the band held before nor on the cells beyond it, nor on the order the cells are visited in. The variables extended are the conserved ones, so nothing in it names a gas.

3. Velocity. Every cell takes the velocity of its owner, from the states of the stage, into one field; the interface is moved with it, so the level set sees the velocity of the flow at the order of the integrator.

4. Operators. Each material applies its own flux operator to its whole field, and the output is set to zero beyond its band, where nothing reads. The time step is the smallest of the materials’, each looking at the cells it owns.

5. Interfaces. Each level-set function is advected with the velocity field by the operator of the level set.

After the step the functions are made to agree with the owner map (the fix-up of Losasso, Shinar, Selle and Fedkiw 2006): in every cell, with $a$ the largest function and $b$ the second largest,

\phi_m\;\leftarrow\;\phi_m-\tfrac12\,(a+b),

so that the owner’s is $(a-b)/2\ge0$ and every other is at most $(b-a)/2\le0$; an overlap and a gap are repaired by the same shift, the interface lies where two functions are equal, the differences between functions are kept, and two functions that are opposites are not changed. Applying it twice changes no owner and no sign and the values by a few units in the last place at most.

What is checked

One material is Solver

The time step and the checks on a step are the code of Solver, shared (crates/merlin-finite-volume/src/solver/limits.rs). One material in the coupled solver gives the bits of Solver, after each of the first steps and at the end, with the checkpoint contents, at 1, 2, 4 and 8 threads, for Sod’s problem with each integrator, a fast inflow that the time step must see, a 2D blast and a 3D flow with walls (unsplit and Strang), a cylindrical blast with an axis, and Lax-Wendroff with Strang splitting; the refusals are Solver’s.

Enforced by crates/merlin-multimaterial/tests/one_material/identical.rs: sod_in_one_dimension_with_each_integrator, inflow_faster_than_the_domain_limits_the_step_as_it_does_in_solver, a_cylindrical_blast_with_radial_weights_and_an_axis, a_blast_in_two_dimensions_unsplit_and_strang, a_smooth_flow_in_three_dimensions_unsplit_and_strang, the_comparison_can_fail; crates/merlin-multimaterial/tests/one_material/lax_wendroff.rs: with_strang_and_forward_euler_it_is_the_run_of_solver, it_is_refused_where_solver_refuses_it; crates/merlin-multimaterial/tests/one_material/errors.rs: the_refusals_of_solver_are_the_refusals_of_the_coupled_solver.

Owners, layers and the fix-up

The layers are compared with the exact face-to-face distance to the nearest cell a material owns, for periodic and non-periodic axes and widths 0, 1, 2 and 5; the owner map, the layers and the fix-up are the same bits for four block layouts and 1 to 8 threads; after the fix-up exactly one function is not negative in every cell and it is the owner’s, for three materials with overlaps, gaps and exact ties.

Enforced by crates/merlin-multimaterial/tests/level_sets/layers.rs: a_layer_is_the_number_of_steps_from_a_cell_the_material_owns; crates/merlin-multimaterial/tests/level_sets/owners.rs: the_owner_of_a_cell_is_the_function_that_is_largest_there_and_a_tie_goes_to_the_lower_index, the_owner_map_does_not_depend_on_the_blocks_or_the_threads; crates/merlin-multimaterial/tests/level_sets/fix_up.rs: after_the_fix_up_every_cell_has_exactly_one_function_that_is_not_negative_its_owners, an_overlap_and_a_gap_are_repaired_by_the_same_shift_and_the_interface_sits_where_two_are_equal, applying_the_fix_up_twice_changes_nothing_but_the_last_bit, two_functions_that_are_opposites_are_not_changed_at_all, a_tie_leaves_exactly_one_function_that_is_not_negative, the_fix_up_does_not_depend_on_the_blocks_or_the_threads.

The extension into the band

Across a plane interface the state continues along the normal to $10^{-13}$, the real cells and the cells beyond the band are not changed, and two different garbages in the band and beyond give the same bits. Round a circle, with a state that depends on the angle alone (so that its exact extension is itself), the error falls at first order, as the rule is:

Cells	48	96	192
Largest error in the band	0.0328	0.0213	0.0114
Order of the pair		0.62	0.90

Enforced by crates/merlin-multimaterial/tests/two_materials/extension.rs: across_a_plane_interface_the_state_continues_along_the_normal_exactly, what_the_band_held_before_does_not_matter, the_extension_does_not_depend_on_the_blocks_or_the_threads, round_a_circle_the_extension_converges_to_the_state_that_is_constant_along_the_normal.

Preparing the interface physics (9.3, step 4). The solver now asks the ghost fill to fill every material at once (GhostFill::fill_all, with a FillRequest), once per stage, so that the ghost state of one material can use the real state of the other; the default is the extension above, material by material, to the last bit. A slot can lend its material at an interface point (MaterialSlot::with_material_at_interface): from a real state interpolated to the point and the normal that points out of the material, a gas builds its response in the frame of that normal and the ghost cell it writes from the traction and normal velocity the interface imposes (the pressure $-t^*$, the velocity $u^*$ along the normal, the gas’s own tangential velocity and the density behind its own wave); a scalar law has none, and keeps the extension. Which law, solver and response to separation (Fail between fluids, FreeSurface for solids, which is phase 10) apply to a pair of materials is a table by name (PairLaws). Nothing here changes a result yet.

The star values and the ghost states (9.3, steps 5 and 6). InterfaceStates is the ghost fill that gives each ghost cell its state from the interface problem. At each layer-1 ghost cell of a material (the inside) it finds the foot point on the interface, $x-\phi\nabla\phi/|\nabla\phi|^2$, and the normal $n=-\nabla\phi/|\nabla\phi|$ out of the material (the functions are positive inside); samples the conserved state of each material one cell width from the foot point into its own region, interpolated from its real cells only (multilinear, widened to 1.5 and 2 cell widths where the stencil would reach another material, with a typed error if it still cannot); lends both materials, one inside the other; solves the interface problem of their pair once; and keeps $t^*$, $u^*$ (along the normal out of the inside) and $n$ in scratch fields allocated once. The material across the interface is the owner of the cell at the far probe. The inside material writes its ghost cell from its side of the solution, inside the same lend: the state behind its own wave at $t^*$ and $u^*$ (the pressure $-t^*$, the density from its wave curve, the velocity $u^*n$) plus its own tangential velocity, so that a vortex sheet is not smeared by the other material’s. Every ghost state is computed and kept before any is written (compute, then apply); then layers 2 and beyond are filled outward from layer 1 by the extension of 9.2, so the whole band carries the star state along the normal. A cell with no response to offer keeps the extension of its own state. The probe reads real cells across the boundaries of the blocks and nothing is written while it does, so the star values are the same bits for any blocks and threads, and cells a material does not own may hold anything.

What is checked. For two gases of one kind (Sod, the 123 problem, two shocks, a strong blast) the star values at the ghost cells of both sides are those of the exact Riemann solver to $10^{-12}$, the other side getting the same traction and the opposite velocity; between air and water they are those of the direct solve. Round a circle, with states that vary linearly in space, they are those of the probe states of a geometry written out in the test, to $10^{-10}$: a probe at another distance or a normal of the wrong sign would differ by far more. The ghost states. For the same four Riemann problems, with Exact, the ghost cells of the left material (layers 1 to 4) hold the exact star state on the left of the contact and those of the right material the one on its right, to $10^{-12}$. Air and water at one pressure and one velocity round a disc (every direction of the interface, an impedance ratio of about 4000) keep their own states in every layer with every solver. With different velocities, each ghost cell has the interface’s normal velocity and its own gas’s tangential one. Sod as two materials of one gas, with the interface ahead of the shock, on the contact and in the rarefaction fan, at 100, 200 and 400 cells to $t=0.2$, for Hllc and Exact: the error of the density against the exact solution is at most 1.64 times that of the single-material solver (0.81 to 0.96 on the contact), converges at an observed order of about 0.9 (the single-material solver 0.88), and the difference between the two runs falls by a factor of 2.7 to 4.2 from 100 to 400 cells; the extension of 9.2, which does not transmit the waves, is 5 to 100 times worse than one material. Courant slack: for a gas at $10^9$ Pa against water at $10^5$ Pa, and the other way round, the fastest wave in the layer-1 ghost cells is at most 1.01 times the one that the time step scans in the cells the materials own, against a limit of 2.5 at a Courant number of 0.4, so the scan of the real cells stays. The velocity of the interface. The interface is a material surface, so it moves at the normal velocity $u^*$ that both sides share, not at the velocity of the real cell next to it, which lags when a wave arrives. GhostFill::interface_velocity replaces, in the cells next to an interface, the velocity the level sets are advected with by $u^*n$ plus the tangential velocity of the cell’s owner (a slip interface does not move with the tangential velocity of either side), carried along the normal through the band. For Sod’s problem as two materials with the interface at the diaphragm, the interface ends within 0.02 to 0.16 of a cell of the exact $0.5+u^*t$ at 100, 200 and 400 cells, against 0.84 to 1.2 cells with the velocity of the owner (the control, the same ghost states). What it costs (release build, one thread; examples/interface_cost.rs): building a material’s side at a point is 8 to 9 ns; the solve itself is 22 ns for Acoustic, 91 ns for Hllc and about 810 ns for Exact; the whole stage, with the probe and the sweep of the grid, adds about 1 to 2.5 µs per ghost cell (timings noisy to 30 %) to the extension, which is 15 to 40 % of its time on a disc of 128 to 512 cells a side; the solve is a small part of it. Limits: an interface that runs into the edge of the domain, a corner of a region or a region thinner than the stencil is the error CannotInterpolate, naming the cell; nothing wraps round a periodic boundary.

Enforced by crates/merlin-multimaterial/tests/interface_states/lender.rs: a_gas_slot_lends_its_response_in_the_frame_of_its_outward_normal, the_normal_velocity_changes_sign_with_the_normal, a_ghost_cell_written_from_the_states_own_values_is_the_state, a_ghost_cell_has_the_star_pressure_and_normal_velocity_and_its_own_tangential_velocity, a_scalar_slot_has_no_response_and_an_inadmissible_probe_is_an_error, two_materials_lent_one_inside_the_other_give_the_ghost_cell_of_the_direct_solve; crates/merlin-multimaterial/tests/interface_states/pair_laws.rs: a_pair_is_unordered_and_the_others_keep_the_default; crates/merlin-multimaterial/tests/interface_states/stars.rs: the_star_values_of_one_gas_are_those_of_the_exact_riemann_solver, the_star_values_between_air_and_water_are_those_of_the_direct_solve, star_values_round_a_circle_are_those_of_the_probe_states_of_the_geometry, the_star_values_do_not_depend_on_the_blocks_or_the_threads, what_a_material_holds_where_it_does_not_own_the_cell_changes_nothing, a_scalar_material_next_to_a_gas_has_no_star_values, a_probe_that_cannot_be_filled_is_an_error_that_names_the_cell; crates/merlin-multimaterial/tests/interface_states/ghost_states.rs: the_ghost_cells_hold_the_exact_star_state_on_their_own_side_of_the_contact, a_contact_between_air_and_water_at_any_angle_leaves_the_ghost_states_unchanged, a_ghost_cell_has_the_normal_velocity_of_the_interface_and_the_tangential_velocity_of_its_own_gas; crates/merlin-multimaterial/tests/interface_states/sod.rs: sod_with_the_interface_ahead_of_the_shock_is_as_accurate_as_one_material, sod_with_the_interface_on_the_contact_is_as_accurate_as_one_material, sod_with_the_interface_in_the_rarefaction_is_as_accurate_as_one_material; crates/merlin-multimaterial/tests/interface_states/velocity.rs: the_interface_follows_the_contact_speed_better_than_the_velocity_of_the_owner; crates/merlin-multimaterial/tests/interface_states/courant.rs: the_ghost_states_of_a_strong_shock_stay_within_the_speed_the_time_step_scans; crates/merlin-multimaterial/tests/two_materials/fill_all.rs: the_solver_fills_all_the_materials_in_one_call_per_stage_and_gives_the_same_bits.

Two gases from a problem file

A problem file with several [[material]] tables runs through the coupled solver (phase 9.3, step 10): the regions give the owner of each cell, the level-set function of a material is the signed distance to the cells it owns (sign, then fast sweeping, Zhao 2005), an optional [[coupling]] table picks the solver of a pair (hllc, the default, or exact), and the gases use the scheme of [numerics] with the HLLC flux. Two shipped problems, at 400 cells: 1d/multimaterial/gas_gas.toml (air against a monatomic gas, Sod’s states, $t=0.2$) and 1d/multimaterial/gas_water.toml (water at $10^9$ Pa against air at $10^5$ Pa and 50 kg/m³, $t=2\times10^{-4}$ s). In both the interface ends within half a cell of $0.5+u^*t$ with $u^*$ of the exact interface problem (a hundredth of a cell for the gases, 0.27 of a cell for the strong jump of water into air), and the cell on each side of it holds the star pressure and velocity to 1 %. The same bits come out for one block on one thread and five blocks on four threads, and the HDF5 and XDMF files are written with the variables tagged by material.

Enforced by crates/merlin-cli/tests/run/multimaterial_problems.rs: two_gases_in_a_shock_tube_end_with_the_interface_where_the_exact_solve_puts_it, a_shock_from_water_into_air_ends_with_the_interface_where_the_exact_solve_puts_it, two_gases_give_the_same_bits_for_any_blocks_and_threads, a_two_gas_run_writes_the_hdf5_and_xdmf_files.

Other geometries (step 10b). The same files run in two dimensions, on an axisymmetric (r, z) grid and in spherical symmetry; the geometry belongs to the grid, and the level sets, the probe and the interface problem work in the coordinates of the flat (r, z) plane. A bubble of a monatomic gas in air in the uniform flow (0.7, 0.4) (2d/multimaterial/bubble_in_flow.toml) and a ring in the axial flow (0, 0.5) (cylindrical/multimaterial/ring_in_flow.toml) keep the pressure and the velocity uniform to $10^{-12}$ over every cell either material owns, with no radial velocity for the ring, and are carried to within a quarter of a cell of $u\,t$. A ball of air at 5000 atmospheres in water in spherical symmetry (spherical/multimaterial/gas_ball_in_water.toml, $t=10^{-4}$ s) pushes the interface from $r=0.3$ to 0.3267 and loses $2\times10^{-4}$ of the total mass, the non-conservation of a ghost-fluid method at the interface; the bound in the test is a regression of that measurement, not a comparison with a reference (9.4). In three dimensions the ghost states round a ball of air in water at one pressure and one velocity are the materials’ own in every layer, for every solver. Found on the way: the level sets of the problem file started from the bare signs of the owners, which the fast-sweeping reinitialisation turned into zeros next to the interface; the fix-up then gave those cells to the other material, so each region shrank by a layer, and a bubble sliding along an axis of the grid stopped the run with CannotInterpolate. They now start from the signs averaged over the $3^D$ neighbours, and the owners at the start are exactly the cells of the regions. Limits: an interface within about two cells of the axis of symmetry (or of the centre of a sphere) is CannotInterpolate naming the cell, because the probe does not mirror across it; the interpolation of a probe is multilinear in the coordinates, not weighted by the volume of the cells, which is first-order consistent and not exact on a strongly curved grid; and gases with swirl cannot be in contact (refused by the configuration).

Enforced by crates/merlin-cli/tests/run/multimaterial_geometry.rs: a_bubble_in_a_uniform_flow_keeps_the_pressure_and_the_velocity_and_is_carried, a_bubble_moving_along_an_axis_or_a_diagonal_runs_and_keeps_the_pressure, the_owners_at_the_start_are_the_cells_of_the_regions, a_ring_in_an_axial_flow_keeps_the_balance_of_the_geometric_terms, an_interface_on_the_axis_is_an_error_that_names_the_cell, a_gas_ball_in_water_pushes_the_interface_out_and_keeps_the_mass; crates/merlin-multimaterial/tests/interface_states/ghost_states.rs: a_contact_round_a_ball_in_three_dimensions_leaves_the_ghost_states_unchanged; crates/merlin-cli/src/config/tests.rs: gases_with_swirl_cannot_be_in_contact_yet.

A bubble in uniform flow

Two ideal gases ($\gamma=1.4$ and $5/3$, densities 1 and 0.3) at the same pressure and velocity $(0.7,0.4)$ on a periodic square, separated by a smooth periodic function; the exact solution is the bubble carried unchanged. Every cell a material owns or reads stays at its state to $10^{-12}$, and the function of the bubble follows the exact translation (to $t=0.2$, SSP-RK3, MUSCL with HLLC):

Cells	24	48	96
Largest error of $\phi$	$1.98\times10^{-4}$	$6.23\times10^{-6}$	$1.95\times10^{-7}$
Order of the pair		4.99	5.00

At least second and a half order is required; about fifth is observed because the time error of SSP-RK3 is small for a translation. Strang splitting moves the bubble as well as the unsplit method.

Each gas carries a smooth pattern in its density in the figure, so the real cells, the ring of ghost cells (the extension along the normal) and the hole where nothing reads are all visible; the pressure and the velocity are uniform to $10^{-15}$.

Enforced by crates/merlin-multimaterial/tests/two_materials/bubble.rs: both_gases_stay_as_they_were_and_the_bubble_moves_with_the_flow, the_interface_follows_the_flow_at_the_order_of_the_integrator_or_better, the_split_method_moves_the_bubble_as_well, every_cell_is_moved_by_the_velocity_of_its_owner, the_answer_does_not_depend_on_the_blocks_or_the_threads.

The order in time of the coupling

In one dimension a pressure wave accelerates the interface of two samples of one gas from rest, by a fifth of a cell in the run, so that the owner map does not change and the whole coupled system is a smooth function of time. The error of $\phi$ against a run with a step 64 times smaller falls with a fitted order of 2.91 (SSP-RK3 is third order):

Courant number	0.4	0.2	0.1	0.05
Largest error of $\phi$	$1.17\times10^{-7}$	$1.34\times10^{-8}$	$2.50\times10^{-9}$	$2.49\times10^{-10}$

A coupling that read the velocity once per step would have order one: the lagged oscillator of merlin-core shows it.

Enforced by crates/merlin-multimaterial/tests/two_materials/acoustic.rs: the_coupling_of_the_flow_and_the_interface_has_the_order_of_the_integrator; crates/merlin-core/src/time/stepping.rs: coupled_fields_have_the_order_of_the_integrator, coupling_at_each_stage_beats_coupling_once_per_step.

Cells that nothing reads, and any blocks and threads

With every cell of each state beyond its band set to NaN, eight steps give the same bits in all the cells each material owns and in both functions as the clean run. The whole coupled state is the same bits for four block layouts and 1, 2, 4 and 8 threads, unsplit and Strang.

Enforced by crates/merlin-multimaterial/tests/two_materials/acoustic.rs: the_cells_beyond_the_band_are_never_read; crates/merlin-multimaterial/tests/two_materials/bubble.rs: the_answer_does_not_depend_on_the_blocks_or_the_threads.

Checkpoint and restart

The checkpoint of a coupled run holds the state of each material and each function (file version 2; version 1, one field, is still read). A run stopped part-way, saved to an HDF5 file and continued on a different number of threads ends in the bits of the run that was not stopped: the states, the functions, the owner map, the layers and the clock, unsplit and Strang. Restoring the states without the functions gives a different run, and a checkpoint of other fields is refused with the field named.

Enforced by crates/merlin-multimaterial/tests/two_materials/restart.rs: a_coupled_run_saved_to_a_file_and_continued_ends_in_the_bits_of_one_that_was_not, restoring_the_states_without_the_functions_would_not_give_the_same_bits, a_checkpoint_of_another_set_of_fields_is_refused_and_leaves_the_solver_alone; crates/merlin-io/tests/checkpoint/main.rs: several_fields_read_back_bit_for_bit_in_order, a_version_1_file_is_read_as_a_checkpoint_of_one_field; crates/merlin-driver/tests/drive/coupled_fields.rs: a_coupled_run_stopped_and_restarted_ends_in_the_bits_of_one_that_was_not; crates/merlin-cli/tests/run/restart.rs: a_checkpoint_written_before_version_2_is_still_restarted_exactly.

What a coupled run writes, and what it conserves

For each snapshot: the state of each material in its natural variables (rho@ambient, u@ambient, …) over the whole grid; the function of each material; the owner of each cell; the layer of each cell for each material; and the velocity of the owners. The ghost cells in use are shown and the others hidden: the band of each material is refilled from its real cells as the operator would read it at the next stage, and every cell beyond the band is NaN (a hole in a plot, so a material leaves no trail). The integrals are of the conserved variables of each material over the cells it owns, so the mass is not counted in ghost cells.

The extension is not conservative at the interface, and the report says so instead of hiding it. In the bubble in uniform flow the mass of each material over its real cells changes only as the number of cells it owns does, which depends on how the interface crosses the grid: at 48 by 48 cells and $t=0.2$ by +0.20 % (ambient), -1.21 % (bubble) and +0.14 % in total; at 128 by 128 cells and $t=1$ (the example, with a pattern in the density) by -0.009 % and +0.058 %. These are measured, not asserted to be zero; the tests bound them at 2 % only as a guard against a gross error. The conservation error of a real interface problem is measured with the physics of 9.3 in 9.4.

Enforced by crates/merlin-multimaterial/tests/two_materials/output.rs: the_fields_written_are_the_ones_the_names_say, the_ghost_cells_in_use_are_shown_and_the_cells_nobody_reads_are_hidden, the_velocity_written_is_the_one_of_the_owner_of_every_cell, the_mass_of_each_material_is_the_mass_of_the_cells_it_owns.

Cost

A step of two materials costs more than two single-material steps: at 256 by 256 cells with MUSCL and HLL and SSP-RK3 (crates/merlin-multimaterial/examples/two_materials_timing.rs, median of three) it takes 3.80 times the step of Solver on one gas at one thread (1.890 s against 0.498 s) and 5.95 times at eight (0.914 s against 0.154 s); two gases whose operators run over the whole grid would take twice. The owner map, the layers and the ghost fill each sweep the whole grid for a band that is a thin ring, and a stage has many small parallel regions. It is not profiled; the remedy is planned in the roadmap (9.2b). With one material the coupled solver is as fast as Solver (+0.1 %, -0.4 % and +2.2 % in three runs of crates/merlin-multimaterial/examples/one_material_timing.rs).

References

• Aslam, T. D. (2004), A partial differential equation approach to multidimensional extrapolation, Journal of Computational Physics 193, 349-355.
• Fedkiw, R. P., Aslam, T., Merriman, B. and Osher, S. (1999), A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method), Journal of Computational Physics 152, 457-492.
• Losasso, F., Shinar, T., Selle, A. and Fedkiw, R. (2006), Multiple interacting liquids, ACM Transactions on Graphics 25 (SIGGRAPH), 812-819.
• Shu, C.-W. and Osher, S. (1988), Efficient implementation of essentially non-oscillatory shock-capturing schemes, Journal of Computational Physics 77, 439-471.