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 theuserdef_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).
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. Thebuffer
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 theamr_module.f90
do the memory allocations for you) the further course ofradmc3d
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 theradmc3d
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
orlines_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 theexamples/
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 theradmc3d.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
ordust_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 theuserdef_parse_main_namelist()
routine and/or in theuserdef_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 theuserdef_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 thewrite_grid_file()
subroutine from withinuserdef_writemodel()
, then you do not have to explicitly type./radmc3dwritegridfile
as well. Note also thatradmc3d
will automatically call thewrite_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 fileamr_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 theuserdef_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:
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.