Very different multifluid behavior at different resolutions

[emprice]

I'm simulating a disk with FARGO3D and have run into some behavior I can't explain. Working in cylindrical coordinates, my default (low) resolution is 128 radial cells x 384 azimuthal cells. To check convergence, I've run the same simulation at 2x resolution in both dimensions. I would expect more structure and detail but overall similar behavior, and that's not what I'm seeing.

At low resolution, a Jupiter-mass embedded planet clears a very obvious gap, so much so that the accretion streams onto the planet aren't visible in this plot of log surface density of dust 5 (inverse Stokes number = 50). (Apologies for the plot quality, these were just for debugging. x-axis is azimuthal cell number and y-axis is radial cell number, with log surface density in color.)

At high resolution, the same Jupiter-mass embedded planet does not obviously clear its gap, there is much more structure and material caught at the Lagrange points.

I'd appreciate any insight you can provide. I've attached the run parameters, options file, and boundary conditions files in an archive, in case that's helpful to see what might be the issue. The only ideas I have so far:

• I'm using a much lower value of alpha (1e-5) for this alpha-disk than the default. Could that be causing some instability (not enough diffusion?)
• I'm using the mass taper option to mitigate shocks caused by introducing the planet at its final mass. Some of my experimentation suggests this option makes a difference in the plots above, but I can't think of why it should (the option is only used at the point in the code where the potential is computed and nowhere else)

fargo_multifluid.tar.gz

[pbllambay]

Hi Ellen,

I ran your setup with the suggested parameters and resolution and I could not find anything strange. I did not include the mass taper so I get a stronger potential and more features to compare.

On a fresh clone, my execution lines are:

make SETUP=fargo_multifluid
./fargo3d -o "nx=384 ny=128 masstaper=0 outputdir=outputs/low ntot=100" setups/fargo_multifluid/fargo_multifluid.par
./fargo3d -o "nx=768 ny=256 masstaper=0 outputdir=outputs/high ntot=100" setups/fargo_multifluid/fargo_multifluid.par
and the fields dust5dens50.dat for each resolution are:

where the plot line is:

plt.imshow(log10(rho), origin="lower", vmin=-7, vmax=-4.5)

Do your results agree with my plot? Is it possible that you changed something else?

Best,
Pablo

[emprice]

Hi Pablo! Thanks for checking. Spoke offline to @krappleo who suspects my low value of alpha may be leading to RWI. Other than using alpha = 1e-5 and epsilon = 1e-2, my working directory is basically a clean clone of the latest FARGO3D also. But I'm not even getting the same result as what you posted. This is really strange... When I increase alpha to 1e-3, I get something closer to what you posted above, though still with some "extra" features that I don't see in your output. Any ideas for moving forward? (I really need to figure this out because my paper draft can't be resubmitted until we address this comment from a referee.)

[pbllambay]

1. To confirm that we both are running the same thing, could you please use the lines I posted above to run the code and post the resulting images here? These runs just take a few minutes on a laptop.

2. Regarding the viscosity value, it's not advisable to use low resolution with low physical viscosity. If physical viscosity is too low, solutions obtained at low resolution are likely to be dominated by numerical diffusion (truncation errors) rather than physical diffusion. Consequently, low-resolution results might not be comparable to those obtained when the mesh resolution is high enough to allow physical diffusion to play a significant role. Therefore, higher alpha values can help mitigate resolution-dependent differences by ensuring that physical viscosity dominates the solution over numerical diffusion. In any case, if the differences between results obtained at different resolutions are substantial, the reliability of the results should be questioned.

[emprice]

I may be able to shed more light on this now. Since you're showing timestep 50 in your plots, I did the same, with alpha = 1e-5. After just 5 orbits, my results look similar to yours, though with more waves thrown off the planet than in yours. The plots I showed in the original issue post were made after 2000 orbits. Is there any chance your simulations used a different alpha or epsilon? When I made plots with alpha = 1e-3, they looked a lot more like what you posted above.

[pbllambay]

To avoid comparing different things I wrote a simple script that automates the process of cloning FARGO3D into a 'test' directory, downloading your specified setup, compiling the code, and running two simulations at different resolutions. Finally, it uses matplotlib to plot dust5dens50.dat.

Please check that you can get my results now. After that, you can run the simulations for longer time and compare. In any case, my suggestion would be to test convergence with resolution of any of your results. Going for low viscosity necessarily requires a resolution that is high enough such that numerical diffusion does not dominate the dynamics of the phenomena you are measuring. Note that in the limit of zero viscosity, convergence or steady-state is not guaranteed as fine structure can appear on smaller and smaller scales. In that case, the solution would be dominated by numerical diffusion (and by the methods used to solve the hydrodynamics).

import numpy as np
import matplotlib.pyplot as plt
import os

runtest=True
plottest=True

testdir = "test"
repo = "git@github.com:FARGO3D/fargo3d.git"
reponame = "fargo3d"
setupfile = "https://github.com/user-attachments/files/17062162/fargo_multifluid.tar.gz"
exec_low = 'mpirun -np 4 ./fargo3d -o "nx=384 ny=128 masstaper=0 outputdir=outputs/low ntot=100" setups/fargo_multifluid/fargo_multifluid.par'
exec_high = 'mpirun -np 4 ./fargo3d -o "nx=768 ny=256 masstaper=0 outputdir=outputs/high ntot=100" setups/fargo_multifluid/fargo_multifluid.par'

if runtest:
    os.mkdir(testdir)
    os.chdir(testdir)
    os.system("git clone "+repo+" "+reponame)
    os.chdir(reponame)
    os.system("rm setups/fargo_multifluid/*")
    os.system("wget -P setups/ https://github.com/user-attachments/files/17062162/fargo_multifluid.tar.gz")
    os.system("tar -xf setups/fargo_multifluid.tar.gz --directory setups/fargo_multifluid")
    os.system("make SETUP=fargo_multifluid PARALLEL=1")
    os.system(" SETUP=fargo_multifluid")
    os.system(exec_low)
    os.system(exec_high)
else:
    os.chdir(testdir)
    os.chdir(reponame)

if plottest:
    fig = plt.figure(figsize=(8,8))
    ax1 = fig.add_subplot(211)
    ax2 = fig.add_subplot(212)
    for d,ax in zip(["outputs/low/", "outputs/high/"],[ax1,ax2]):
        x = np.loadtxt(d+"domain_x.dat")[:-1]
        y = np.loadtxt(d+"domain_y.dat")[3:-4]
        nx = len(x)
        ny = len(y)
        rho = np.fromfile(d+"dust5dens50.dat").reshape(ny,nx)
        ax.set_title(d)
        ax.imshow(np.log10(rho),aspect="auto",origin="lower",vmin=-7, vmax=-4.5)
    plt.show()