*************************** Setting Up Your Own Problem *************************** Castro problems are organized loosely into groups describing their intent (e.g., science, hydro tests, ...). These groups are sub-directories under the ``Castro/Exec/`` directory. Each problem is then placed in a sub-directory of the appropriate group (for example, ``Castro/Exec/hydro_tests/Sedov`` holds the Sedov test problem). To create a new problem, you will create a new directory under one of the groups and place in it the following files: * ``GNUmakefile`` : the makefile for this problem. This will tell Castro what options to use and what network and EOS to build. * ``problem_initialize.H`` and ``problem_initialize_state_data.H`` : this holds the problem initialization routines. MHD and radiation problems require an additional file. * ``_prob_params`` (optional) : a list of runtime parameters that you problem will read. These parameters are controlled by the set in the inputs file as ``problem.param`` where ``param`` is one of the parameters listed in ``_prob_params``. * ``Make.package`` : this is a makefile fragment that is included during the build process. It tells the build system about any problem-specific files that need to be compiled. * ``inputs`` : this is the main inputs file that controls Castro and AMReX's behavior. The best way to get started writing a new problem is to copy an existing problem setup into a new directory. .. index:: _prob_params Runtime Parameters ------------------ The problem-specific runtime parameters are defined in ``_prob_params``. This has the form:: name datatype default namelist? size Here: * `name` is the name of the variable to put into ``probdata_module`` * `datatype` is one of ``real``, ``integer``, ``character``, or ``logical``. * `default` is the default value of the runtime parameter. It may be overridden at runtime by reading from the namelist. * `namelist` indicates if the variable should be able to be set as a runtime parameter. If it is empty or marked with ``n``, then the variable will still be creates in the ``problem`` namespace, but it will not be able to be set via the commandline or inputs file. A common usage of this is to define global variables that might be set at problem initialization that are used elsewhere in Castro. * `size` is for arrays, and gives their size. It can be any integer or variable that is known to the ``probdata_module``. If you need a variable from a different module to specify the size (for example, ``nspec`` from ``network``), then you express the size as a tuple: ``(nspec, network)`` and the appropriate use statement will be added. The variables will all be initialized for the GPU as well. Problem Initialization ---------------------- Here we describe the main problem initialization routines. .. index:: initialize_problem * ``initialize_problem()`` This C++ routine is primarily responsible for doing any one-time initialization for the problem (like reading in an initial model through the ``model_parser.H`` functionality—see the ``toy_convect`` problem setup for an example. This routine can also postprocess any of the parameters defined in the ``_prob_params`` file, which are defined in ``problem`` namespace. .. note:: many problems set the value of the ``problem::center[]`` array from the ``problem`` namespace. This is used to note the center of the problem (which does not necessarily need to be the center of the domain, e.g., for axisymmetric problems). ``center`` is used in source terms (including rotation and gravity) and in computing some of the derived variables (like angular momentum). If you need coordinate information, it can be obtained by constructing a ``Geometry`` object using ``DefaultGeometry()`` and accessing its ``ProbLo()`` and ``ProbHi()`` methods. * ``problem_initialize_state_data()``: This routine will initialize the state data in a given zone. The arguments passed are: - ``i``, ``j``, ``k``: the index of the zone to fill the data in - ``state``: an array containing the simulation state data - ``GeomData``: a ``GeometryData`` object that can be used for obtaining ``dx``, ``problo``, ``probhi``, etc. Filling data is done by simply writing to ``state(i,j,k,URHO)``, etc. .. _create:bcs: Boundary Conditions ------------------- .. index:: boundary conditions Standard boundary conditions, including outflow (zero-gradient), periodic, and symmetry (reflect) are handled by AMReX directly. Castro has a special hydrostatic boundary condition that can be used for the lower boundary. It is accessed by setting the ``castro.lo_bc`` flag to 1 in the vertical coordinate direction, e.g., for 2-d as:: castro.lo_bc = 0 1 The flag value 1 is traditionally named "inflow" by AMReX, but generally means that the boundary implementation is left to the user. To tell Castro to use the hydrostatic boundary condition here, we set:: castro.yl_ext_bc_type = 1 castro.hse_interp_temp = 1 castro.hse_reflect_vels = 1 The first parameter tells Castro to use the HSE boundary condition for the lower y direction. In filling the ghost cells, hydrostatic equilibrium will be integrated from the last interior zone into the boundary. We need one more equation for this integration, so we either interpolate the density or temperature into the ghost cells, depending on the value of ``castro.hse_interp_temp``. Finally, ``castro.hse_reflect_vels`` determines how we treat the velocity. The default is to give is a zero gradient, but in tests we've found that reflecting the velocity while integrating the HSE profile can be better. For modeling a plane-parallel hydrostatic atmosphere, using the hydrostatic boundary conditions instead of a simple symmetry boundary is essential when using the standard CTU PPM solver. A different special boundary condition, based on outflow, is available at the upper boundary. This works together with the ``model_parser`` module to fill the ghost cells at the upper boundary with the initial model data. You set this as:: castro.hi_bc = 2 2 castro.fill_ambient_bc = 1 castro.ambient_fill_dir = 1 castro.ambient_outflow_vel = 1 where ``ambient_fill_dir`` is the 0-based direction to fill using an ambient state defined by the problem setup. In this example, we will override the outflow (2) boundary condition in the y-direction. That problem setup needs to fill the ``ambient_state[:]`` array defined in ``ambient.H``. An example of using this boundary is in the ``flame_wave`` problem. The implementations of these boundary conditions is found in ``Castro/Source/problems/Castro_bc_fill_nd.cpp``. Optional Files -------------- The follow problem-specific files are optional. There are stubs for each of these in the main source tree. - ``problem_checkpoint.H``, ``problem_restart.H`` : These provides two routines, respectively ``problem_checkpoint`` and ``problem_restart`` that can be used to add information to the checkpoint files and read it in upon restart. This is useful for some global problem-specific quantities. For instance, the ``wdmerger`` problem uses this to store center of mass position and velocity information in the checkpoint files that are used for runtime diagnostics. The name of the checkpoint directory is passed in as an argument. - ``problem_tagging.H`` This implements problem-specific tagging for refinement, through a the function ``problem_tagging``. The full hydrodynamic state (State_Type) is passed in, and the problem can mark zones for refinement by setting the tag variable for a zone to set. An example is provided by the ``toy_convect`` problem which refines a rectangular region (fuel layer) based on a density parameter and the H mass fraction. .. _problem_derives: - ``Problem_Derives.H``, ``Problem_Derive.H``, and ``Problem_Derive.cpp`` Together, these provide a mechanism to create derived quantities that can be stored in the plotfile. ``Problem_Derives.H`` provides the C++ code that defines these new plot variables. It does this by adding them to the ``derive_lst``—a list of derived variables that Castro knows about. When adding new variables, a descriptive name, a C++ routine that does the deriving, and component of ``StateData`` are specified. The other two files provide the header and implementation of the function that computes the derived variable. A example is provided by the ``reacting_bubble`` problem, which derives several new quantities (perturbations against a background one-dimensional model, in this case). - ``Prob.cpp``, ``Problem.H`` These files provide problem-specific routines for computing global diagnostic information through the sum_integrated_quantities functionality that is part of the ``Castro`` class. An example is provided by ``toy_flame``, where an estimate of the flame speed is computed by integrating the mass of fuel on the grid. Model Parser ------------ Many problem setups begin with a 1-d initial model that is mapped onto the grid. The ``model_parser.H`` provides the functions that read in the initial model and map it on the Castro grid. To enable this, add:: USE_MODEL_PARSER = TRUE to the problem ``GNUmakefile``. There are 2 other parameters that can be set in the makefile to control the initial model storage: * ``MAX_NPTS_MODEL``: is the maximum number of data points in the 1-d initial model. This needs to be known at compile time so we can make the data managed for GPUs. * ``NUM_MODELS``: this is the number of different initial models we want to managed. Typically we only want 1, but some problems, like ``flame_wave`` use 2, applied to different portions of the domain. The general form of the initial model is:: # npts = 896 # num of variables = 6 # density # temperature # pressure # carbon-12 # oxygen-16 # magnesium-24 195312.5000 5437711139. 8805500.952 .4695704813E+28 0.3 0.7 0 585937.5000 5410152416. 8816689.836 0.4663923963E+28 0.3 0.7 0 The first line gives the number of points in the initial model, the next gives the number of variables, followed by the variable names (one per line), and then the data. The data begins with the coordinate and then the variables in the model, with one data point per line. When the model is read, the variables listed in the file are matched to the ones that Castro knows about. If the variable is recognized, then it is stored in the model data, otherwise, it is ignored. The data can then be mapped onto the grid using the ``interpolate()`` function, e.g., :: Real dens = interpolate(height, model::idens); This fills ``dens`` with the density at the position ``height``. In addition to density, you can specify temperature (``model::itemp``), pressure (``model::ipres``), species (indexed from ``model::ispec``), or an auxiliary quantity (indexed from ``model::iaux``).