Modifying RADMC-3D: Internal setup and user-specified radiative processes

It has been mentioned several times before that as an alternative to the standard compile once-and-for-all’ philosophy, one can also use RADMC-3D by modifying the code directly so that ``radmc3d` will have new functionality that might be of use for you. We refer to Section Making special-purpose modified versions of RADMC-3D (optional) for an in-depth description of how to modify the code in a way that is non-invasive to the main code. We urge the reader to read Section Making special-purpose modified versions of RADMC-3D (optional) first before continuing to read this chapter. In all of the following we assume that the editings to the fortran files are done in the local way described in Section Making special-purpose modified versions of RADMC-3D (optional) so that the original source files in the src/ directory stay unaffected, and only local copies are edited.

Setting up a model inside of RADMC-3D

The most common reason for editing the code itself is for setting up the model internally rather than reading in all data via input files. For a list of advantages and disadvantages of setting models up internally as opposed to the standard way, see Section Some caveats and advantages of internal model setup below. Setting up a model within RADMC-3D is done by making a local copy of the file userdef_module.f90 and editing it (see Section Making special-purpose modified versions of RADMC-3D (optional)). This file contains a set of standard subroutines that are called by the main program at special points in the code. Each subroutine has a special purpose which will be described below. By keeping a subroutine empty, nothing is done. By filling it with your own code lines, you can set up the density, temperature or whatever needs to be set up for the model. In addition to this you can do the following as well:

  • Add new variables or arrays in the module header (above the contains command), which you can use in the subroutines of the userdef_module.f90 module. You are completely free to add any new variables you like. A small tip: it may be useful (though not required) to start all their names with e.g. userdef_ to make sure that no name conflicts with other variables in the code happen.

  • Add new subroutines at will (below the contains command) which you can call from within the standard subroutines.

  • Introduce your own radmc3d command-line options (see Section The pre-defined subroutines of the userdef_module.f90).

  • Introduce your own radmc3d.inp namelist variables (see Section The pre-defined subroutines of the userdef_module.f90).

_images/dataflow-basic-userdef.svg

Fig. 27 Pictographic representation of the dataflow for the case when you define your model internally using the userdef_module.f90.

Often you still want some of the input data to be still read in in the usual way, using input files. For instance, you may want to still read the dustopac.inp and the opacities using the dustkappa_xxx.inp files. This is all possible. Typically, you simply keep the files you still want RADMC-3D to read, and omit the files that contain data that you allocate and set in the userdef_module.f90. This is all a bit complicated, so the best way to learn how to do this is to start from the example directories in which a model is set up with the userdef_module.f90 method.

In Fig. Fig. 27 the dataflow for the user-defined model setup is graphically depicted.

The pre-defined subroutines of the userdef_module.f90

The idea of the userdef_module.f90 is that it contains a number of standard pre-defined subroutines that are called from the main.f90 code (and only from there). Just browse through the main.f90 file and search for the sequence calluserdef_ and you will find all the points where these standard routines are called. It means that at these points you as the user have influence on the process of model setup. Here is the list of standard routines and how they are used. They are ordered roughly in chronological order in which they are called.

  • userdef_defaults()

    This subroutine allows you to set the default value of any new parameters you may have introduced. If neither on the command line nor in the radmc3d.inp file the values of these parameters are set, then they will simply retain this default value.

  • userdef_commandline(buffer,numarg,iarg,fromstdi,gotit)

    This subroutine allows you to add your own command-line options for radmc3d. The routine has a series of standard arguments which you are not allowed to change. The buffer is a string containing the current command line option that is parsed. You will check here if it is an option of your module, and if yes, activate it. An example is listed in the code. You an also require a second argument, for which also an example is listed in the original code.

  • userdef_commandline_postprocessing()

    After the command line options have been read, it can be useful to check if the user has not asked for conflicting things. Here you can do such checks.

  • userdef_parse_main_namelist()

    Here you can add your own namelist parameters that read from the radmc3d.inp file. An example is provided in the original code.

  • userdef_main_namelist_postprocessing()

    Also here, after the entire radmc3d.inp file has been read and interpreted, you can do some consistency checks and postprocessing here.

  • userdef_prep_model()

    This routine can be used if you wish to set up the grid not from input files but internally. You will have to know how to deal with the amr_module.f90 module. You can also set your own global frequency grid here. And finally, you can set your own stellar sources here. In all cases, if you set these things here (which requires you to make the proper memory allocations, or in case of the gridding, let the amr_module.f90 do the memory allocations for you) the further course of radmc3d will skip any of its own settings (it will simply detect if these arrays are allocated already, and if yes, it will simply not read or allocate them anymore).

  • userdef_setup_model()

    This is the place where you can actually make your own model setup. By the time this subroutine is called, all your parameters have been read in, as well as all of the other parameters from the original radmc3d code. So you can now set up the dust density, or the gas velocity or you name it. For all of these things you will have to allocate the arrays youself (!!!). Once you did this, the rest of the radmc3d code won’t read those data anymore, because it detects that the corresponding arrays have already been allocated (by you). This allows you to completely circumvent the reading of any of the following files by making these data yourself here at this location:

    • amr_grid.inp or in the future the input files for any of the other griding types.

    • dust_density.inp

    • dust_temperature.dat

    • gas_density.inp

    • gas_temperature.inp

    • gas_velocity.inp

    • microturbulence.inp

    • levelpop_XXX.dat

    • numberdens_XXX.inp

    To learn how to set up a model in this way, we refer you for now to the ioput_module.f90 or lines_module.f90 and search for the above file names to see how the arrays are allocated and how the data are inserted. I apologise for not explaining this in more detail at this point. But examples are or will be given in the examples/ directory.

  • userdef_dostuff()

    This routine will be called by the main routine to allow you to do any kind of calculation after the main calculation (for instance after the monte carlo simulation). This is done within the execution-loop.

  • userdef_compute_levelpop()

    This is a subroutine that can be called by the camera module for on-the-fly calculation of level populations according to your own recipe. This may be a bit tricky to use, but I hope to be able to provide some example(s) in the near future.

  • userdef_srcalp()

    This subroutine allows you to add any emission/absorption process you want, even fake ones. For instance, you could use this to create nicely volume-rendered images of your 3-D models with fake opacities, which are chosen to make the image look nice and/or insight-giving. You can also use this to add physical processes that are not yet implemented in RADMC-3D. This subroutine allows you full freedom and flexibility to add emissivity and extinction whereever/however you like. To activate it you must set incl_userdef_srcalp=1 in the radmc3d.inp file.

  • userdef_writemodel()

    This allows the user to dump any stuff to file that the user computed in this module. You can also use this routine to write out files that would have been used normally as input file (like amr_grid.inp or dust_density.inp) so that the Python routines can read them if they need. In particular the grid information may be needed by these external analysis tools. Here is a list of standard subroutines you can call for writing such files:

    • write_grid_file()

    • write_dust_density()

    • …more to come…

For now this is it, more routines will be included in the future.

Note that the userdef_compute_levelpop() subroutine, in contrast to all the others, is called not from the main.f90 program but from the camera_module.f90 module. This is why the camera module is the only module that is higher in compilation ranking than the userdef module (i.e. the userdef module will be compiled before the camera module). For this reason the userdef module has no access to the variables of the camera module. For the rest, the userdef module has access to the variables in all other modules.

Note also that not all input data is meant to be generated in this way. The following types of data are still supposed to be read from file:

  • Dust opacity data

  • Molecular fundamental data

Please have a look in the examples/ directory for models which are set up in this internal way.

Some caveats and advantages of internal model setup

Setting up the models internally has several advantages as well as disadvantages compared to the standard way of feeding the models into radmc3d via files. The advantages are, among others:

  • You can modify the model parameters in radmc3d.inp and/or in the command line options (depending on how you allow the user to set these parameters, i.e. in the userdef_parse_main_namelist() routine and/or in the userdef_commandline() routine. You then do not need to run Python anymore (except for setting up the basic files; see examples). Some advantages of this:

    • It allows you, for instance, to create a version of the radmc3d code that acts as if it is a special-purpose model. You can specify model parameters on the command line (rather than going through the cumbersome Python stuff).

    • It is faster: even a large model is built up quickly and does not require a long read from large input files.

  • You can make use of the AMR module routines such as the amr_branch_refine() routine, so you can adaptively refine the grid while you are setting up the model.

Some of the disadvantages are:

  • The model needs to be explicitly written out to file and read into Python or any other data plotting package before you can analyze the density structure to test if you’ve done it right. You can explicitly ask ./radmc3d to call the userdef_writemodel() subroutine (which is supposed to be writing out all essential data; but that is the user’s responsibility) by typing ./radmc3dwritemodel.

  • Same is true for the grid, and this is potentially even more dangerous if not done. You can explicitly ask ./radmc3d to write out the grid file by typing ./radmc3dwritegridfile. Note that if you call the write_grid_file() subroutine from within userdef_writemodel(), then you do not have to explicitly type ./radmc3dwritegridfile as well. Note also that radmc3d will automatically call the write_grid_file() subroutine when it writes the results of the thermal Monte Carlo computation, if it has its grid from inside (i.e. it has not read the grid from the file amr_grid.inp.

  • It requires a bit more knowledge of the internal workings of the radmc3d code, as you will need to directly insert code lines in the userdef_module.f90 file.

Using the userdef module to compute integrals of \(J_\nu\)

With the monochromatic Monte Carlo computation (see Section Special-purpose feature: Computing the local radiation field) we can calculate the mean intensity \(J_\nu\) at every location in the model at a user-defined set of wavelengths. However, as mentioned before, for large models and large numbers of wavelengths this could easily lead to a data volume that is larger than what the computer can handle. Since typically the main motivation for computing \(J_\nu\) is to compute some integral of the the form:

\[Q = \int_0^{\infty} J_\nu K_\nu d\nu\]

where \(K_\nu\) is some cross section function or so, it may not be necessary to store the entire function \(J\) as a function of \(nu\). Instead we would then only by interested in the result of this integral at each spatial location.

So it would be useful to allow the user to do this computation internally. We should start by initializing \(Q(x,y,z)=0\) (or \(Q(r,\theta,\phi)=0\) if you use spherical coordinates). Then we call the monochromatic Monte Carlo routine for the first wavelength we want to include, and multiply the resulting mean intensities with an appropriate \(\Delta\nu\) and add this to \(Q(x,y,z)\). Then we do the monochromatic Monte Carlo for the next wavelength and again add to \(Q\) everywhere. We repeat this until our integral (at every spatial location on the grid) is finished, and we are done. This saves a huge amount of memory.

Since this is somewhat hard to explain in this PDF document, we refer to the example model run_example_jnu_integral/.

STILL IN PROGRESS.

Some tips and tricks for programming user-defined subroutines

Apart from the standard subroutines that must be present in the userdef_module.f90 file (see Section The pre-defined subroutines of the userdef_module.f90), you are free to add any subroutines or functions that you want, which you can call from within the predefined subroutines of Section The pre-defined subroutines of the userdef_module.f90. You are completely free to expand this module as you wish. You can add your own variables, your own arrays, allocate arrays, etc.

Sometimes you may need to know ‘where you are’ in the grid. For instance, the subroutine userdef_compute_levelpop() is called with an argument index. This is the index of the current cell from within which the subroutine has been called. You can now address, for instance, the dust temperature at this location:

temp = dusttemp(1,index)

(for the case of a single dust species). You may also want to know the coordinates of the center of the cell. For this, you must first get a pointer to the AMR-tree structure of this cell. The pointer b is declared as

type(amr_branch), pointer :: b

Then you can point the pointer to that cell structure

b => amr_index_to_leaf(index)%link

And now you can get the x,y,z-coordinates of the center of the cell:

xc = amr_finegrid_xc(b%ixyzf(1),1,b%level)
yc = amr_finegrid_xc(b%ixyzf(2),2,b%level)
zc = amr_finegrid_xc(b%ixyzf(3),3,b%level)

Or the left and right cell walls:

xi_l = amr_finegrid_xi(b%ixyzf(1),1,b%level)
yi_l = amr_finegrid_xi(b%ixyzf(2),2,b%level)
zi_l = amr_finegrid_xi(b%ixyzf(3),3,b%level)
xi_r = amr_finegrid_xi(b%ixyzf(1)+1,1,b%level)
yi_r = amr_finegrid_xi(b%ixyzf(2)+1,2,b%level)
zi_r = amr_finegrid_xi(b%ixyzf(3)+1,3,b%level)

Creating your own emission and absorption processes

RADMC-3D Allows you to add your own physics to the ray-tracing images and spectra. At every point during the ray-tracing process, when it computes the emissivity and extinction coefficients \(j_\nu\) and \(\alpha_\nu\) it calls the userdef_srcalp() subroutine, giving it the index in which cell we are, the frequencies of the different image channels and the src and alp arrays which are for resp.\(j_\nu\) and \(\alpha_\nu\). You can add any process by

src(:) = src(:) + .....
alp(:) = alp(:) + .....

where …… is your formula. You can find the local variables like density and temperature using the index, e.g.:

rho_g = gasdens(index)

You can be completely free in your choices. If you need some information that is not usually read into RADMC-3D, you can add read commands in the userdef_setup_model() subroutine, e.g.:

call read_gas_density(1)

See the example directory examples/run_simple_userdefsrc for more ideas.