Comparison of hydrocode performance in disk-planet interaction

Motivation:

Numerical computation of disk-planet interaction is an increasingly important way to figure out how how the solar and extrasolar planets form. (An easy estimate will convince you that ~1e3 new planetary systems come into being every second, and that's just in the potentially observable Universe, within our event horizon).
It addresses a set of very fundamental processes such as: planet growth, migration in disks, and eccentricity evolution.

This field is rich in opportunities for application of modern CFD techniques. At the same time, experience shows that there is hardly any single CFD method that would be trouble-free under *all* possible circumstances in disk modeling. The list of usual concerns includes but is not restricted to these popular items [send your own list to be included]:

• Boundary conditions influencing the results in a clearly non-physical way, for instance waves bouncing off the solid-wall boundaries and bringing back the angular momentum radiated away from a planet, or causing excessive loss of gas (open boundaries)
• Nonaxisymmetric instabilities of flow of unknown origin at late time (> 100 P, i.e. 100 orbital periods); could be due to physical, growing modes artificially introduced by wave-bouncing boundaries, or purely numerical due to errors connected with the implementations of "fictitious forces" like the Coriolis force in the curvelinear, rotating frames.
• Strange small scale filamentary structure in high-resolution methods run at  high resolution (e.g., PPM run on multi-hundred^2 grids, say, 600^2 or so) in large areas of the disk, resembling some hydrodynamical instability of the disk itself or (in other view) encouraged by the non-stationary flow of gas near the planet (which would then propagate into the whole disk together with the spiral wakes/waves) or (in yet another view!) caused in some way by the algorithm used.
• Code crashes if large planets are introduced instantaneously into disks -- probably due to inability of codes to handle strong shocks in the transient states (some codes, e.g., are not programed to take a step back and repeat the time step which results in violation of the Courant criterion, they merely recommend a shorter time step for the next step. Other difficulties with modeling Jupiter-class bodies may be present as well).
• Planetary wakes or only one of them (outer) often are time-variable. Is this physical or numerical?
• If the flow is stable, it may still be difficult to show that it is an acceptable solution to Navier-Stokes or whatever  form of fluid equations are being solved. Numerical convergence in the limit of infinitely fine mesh may be slow.
• What's more disturbing, and clearly important for the final result of disk-planet interaction, while there may be agreement between the same type of codes, *substantially different* codes might not converge to the same flow. In that context, currently codes seems to group into broadly defined categories of 1. SPH particle hydrodynamiscs, 2. ZEUS-type Eulerian mesh code algorithms, and 3. higher-order schemes like PPM-type and CLAW-type solvers.
• While there is a widespread feeling that some codes provide higher resolution than others, we don't know how to express that quantitatively. How many SPH particles do we need to get the spiral waves resolved to a comparable degree as on a grid of given dimension? Can we compensate for supposedly lower resolution in ZEUS than in PPM (which does more thinking about "subgrid" propagation of disturbances) by simply using a finer mesh, and if so how much finer? Or are some features produced by PPM artifacts caused by some post-shock oscillations or other details?

It would be very valuable to try all of these methods on a unique problem or a set of test problems. It would be even better to have a realistic test problem with analytically known solutions, but we seem to be lacking those in the multidimansional disk-planet interaction.

A number of us, meaning researchers actively doing planet+disk modeling, became interested in conducting the test on a specific disk+planet problem. Garrelt Mellema (Leiden) and Pawel Artymowicz (Stockholm) were considering the idea of a mini workshop in 2003 to do some tests beforehand and compare results during a very specialized meeting. This has been done before in a Santa Barbara program by cosmologists. They jointly published their comparisons or 10 or so codes. Garrelt knows more examples in the field of planetary nebulae modeling.

One practical problem in 2003 is that there are already a number of meetings on planetary systems, so the potential participants may be too busy to attend such a mini-workshop.

The idea was picked up and modified by the steering committee of the new European network on planet formation in Jan. 2003, where Willy Kley (Tubingen), Richard Nelson (London) and Pawel Artymowicz made an effort to formulate a simple setup to be simulated by anybody and everybody.

The RTN Network (coordinated by Andi Burkert, Heidelberg) may spend a little of its networking money on facilitating a meeting, if needed, or in other ways (also, spreading the word about the initiative). But as a minimum, we'll do the work in our spare time separately, and then discuss results via email, convene at meetings such as the Paris meeting in the Summer 2003, and others.

Organizationally, and to facilitate comparison, we will try to process the output data from different codes in one place, which was preliminarily agreed to be Stockholm (other bids are still accepted). The Central  will get ASCII data files,
and produce comparison in accordance with suggestions from participants, maybe even using their proposed software
like IDL scripts, for example.

Here it is, described below, our setup.
Please take a look and comment. Not all the details have been set in stone, there is still
a possibility that we need to change values of parameters etc.
Remember, this isn't exactly a calculation to show us the truth about how Jupiter formed. It's a test of our codes, and
a look into the nature of Keplerian flows perturbed by "planets". We'll think about the direction in which we'll proceed after the first comparisons are done. It would be wonderful to publish if we think the results warrant it.

Last updated: 12 Feb 2003

Will start looking into results:  as soon as available. The idea is to have most data ready for comparison in April 2003.
 

• INPUT/INITIAL CONDITIONS:
• ========================
• 2-D disk  (ok, if you only do 3-d SPH, go ahead, but we need 2-d comparison first),
• with initial flat surface density profile Sigma=const. such that  pi*a^2*Sigma /M_star = 0.002 (if you ever need that; we think you shouldn't initially because we won't initially include self-gravity). a is the semi-major axis of our planet or planetary core, to be precise, and is the unit of distance.
• a=const during the simulation:  a fixed circular orbit, no migration
• willingness to migrate will be expressed as da/dt, computed on the basis of directly measured gravitational torque of the disk on the planet, using 1st order perturbation approach, essentially da/dt/a = 2 * torque/ang. momentum of the core.
• planet's mass = constant, mass ratio  = planet/(planet+star) = mu= {0.001, 0.0001}   (two choices, in order of preference). As some codes may crash with Jupiter, and "Neptune" case of 1/10 Jupiter mass is also very intersesting on its own, we'd like to encourage these two mass parameters.
• locally isothermal equation of state, simulating H/R=0.05 throughout. In other words, soundspeed will be 0.05 times the local Keplerian speed.
• coordinate system - any; rotating or non-rotating frame. Whatever you use, but as many as possible
• calculations on polar grid, first and foremost do a test on uniform azimuthal grid, you may also submit non-uniform grid results. Center grid on the star.
• radial grid: from 0.4a to 2.5a,
• attempt to make square cells near planet.
• initial disk rotation: Keplerian around the star of mass 1.-mu
• ramp up the mass of the planet from 0 to mu smoothly over the first 5 P (orbital periods), using sin(t/5P*pi/2.)^2 ramp.
• apply a smoothed gravitational potential for the planet: -mu/sqrt(r^2 + eps^2) where eps = 0.2 r_L (Roche lobe radius).
• special treatment of boundaries: we'll attempt to eliminate wave reflection on our restricted grid by a simple yet usually effective device. We will interpolate the Sigma(x,y) in the outer and inner disk annulae (0.4-0.5 ad 2.1-2.5)*a as follows: there will be a spatial ramp R(r) decreasing smoothly from 1 at the boundary to 0 on the disk side,
• surface density will be reset after every time step according to the following equation : dX/dt = -(X-X0)/tau*R(r), where X0 is the initial profile of  gas flow quantity (Sigma, and velocity components).
• We don't know exactly at the moment whether to recommend a linear or other R(r); so we recommend linear to the power of 2 (there will be a very smooth transition from disk to the wave-killing zone).
• for comparison, a solid-wall boundary would be good to have, at least a few models.
• in SPH-type codes, if it sounds like a difficult boundary condition to implement, we recommend going to a larger range of r, and forget about the special treetment, we'll see if we can spot any reflected waves or different global evol. of the disk.
• viscosity (artificial or Navier-Stokes): none (if possible), minimum acceptable otherwise (describe)
• we'd like to start with recommendation tau= 1 P(r_BC), or one period of the disk material at the boundary. We'll adjust if it doesn't work. Whoever manages to implement that devise first please report. The idea is to kill the waves but gently, in order to avoid reflections off a hard boundry.
• Those models which use Cartesian grid and a star on it, please smooth the star potential locally, not globally,ie. not through 1/sqrt(r^2 + eps^2) on the whole grid.
• we want to jump by factors of 2 in linear resolution, and continue to as fine meshes as we can, certainly we'll need resolution on the order 0.016666*a near the planet so we propose to depart from that in both directions. But since codes have different diffusivity, that's a weak recommendation, we'll try to compare whatever models you send.
•  
•  
• OUTPUT:
• ========
•  
• at times: t = (2, 5, 10, 20, 50, 100, 200, 500P, ....as far as if you want and manage; suggest to 100P is sufficient)
• ASCII files with a short header, described below,  followed by N lines of data (N=number of cells) containing coordinates and average quantities in the centers of cells: x, y, Sigma, vx, vy   (Cartesian) or, in polar coord:  r, phi, Sigma, vr, vphi.  For SPH, of course, list positions, V's and smoothing lengths of particles
• header should state:  the coord. system, something about the model type, time step used, value of N to simplify reading, and the value of da/dt computed by you (we'll do our estimate to see as well, based on the information about the smoothing of potential that you used)
• please give enough significant digits for reliable computation of derivatives, say, 6 if you compute in single and 12 if in double prec.
• compress the data and provide a link to them on your web page (preferred!) or send them as attachment to the address pawel@astro.su.se. Name the file something like:   Richards_Nirvana.t=2.dat.gz, etc.
• report any crashes, problems.
• remember to indicate to us which aspects of flow you are most interested in. (We'll try to look at Sigma(r) - axisymmetrized density profiles, 2-D density maps, and potential vorticity maps.
• If you do movies, that would be very valuable addition to the output since time-dependence of flow is definitely one of the areas where we see differences and would like to know the reason. Again, we'd encourage posting them on the web for all of us to see.
•  
• All your output will be public domain and can be supplied to all participants on request. (do you want your files linked to this page or not?)
• However, a precondition for participation in this joint program is that we agree not to misuse or maliciously ignore each other's data. (Don't ask me how we're gonna enforce that :-)

Comments welcome,
cheers

EU network "Planets" page
our planet formation group at Stockholm Observatory, Stk Univ.
a little comparison of CLAW and PPM, but not on disk+planet problem