Tutorial: How to add a new parameter¶
- Tutorial: How to add a new parameter
- Add the new parameter into asteroid.inp
- Add the new variable to the appropriate module
- Enable iSALE to read in the new variable
- Use new variable in the source code
- Add parameter information to the database
- Build the code; verify the implementation and documentation
- Add simple test to check whether the input variable is read and used correctly
- Run the tests to verify that changes are good
- Commit changes to repository
This tutorial shows how you can easily extend the functionality of
iSALE. We will show you how to define and insert a new input-parameter, how to update the parameter database, how to use the new parameter in the code and how to write small tests.
iSALE assumes an initially intact material for every layer and the projectile. Here we want to implement the option of defining the default-damage value for each layer and projectile. This can be easily done in quite a few steps as shown below.
Add the new parameter into
For this demonstration, we are using the
Chicxulub-example as provided with
iSALE in the examples-directory. This setup shows the formation of the Chicxulub impact structure in a three-layer target: A thin sedimentary (=_calcite_) target above a granite layer and mantle material composed of dunite. For our new scenario we want to start with a granite layer that is initially partially damaged (damage = 0.5) and a totally damaged sedimentary layer (damage = 1.). This functionality is not yet provided with
iSALE, so we need to implement some changes.
important note: Since every
make installwill overwrite your changes made in the Chicxulub-example directory, it is recommended to prepare the directory structure first. Please read here how this works.
We start with adding some new lines into the
asteroid.inp defining the initial damage of all layers and projectiles. We add the following line to the "
Target Parameters" section of
LAYDAM initial layer damage : 0.D0 : 0.5D0 : 1.D0
We might also want to change the initial damage-state of the projectile, so we add the following line into the "
Projectile ("Object") Parameters" section of
OBJDAM initial object damage : 0.D0
Add the new variable to the appropriate module¶
Now we need to include some modifications in the
iSALE-sourcecode. Before we are able to read in the new variables, we need to define these variables in appropriate modules so that
iSALE knows this variable. The best place for including these "object"-related parameter is in the "object"-structure of
iSALE. This is found in
In this file, you find a so-called "derived datatype"
type type_obj character(len=10) :: name ! material name integer :: mat ! material ID integer :: cppr(0:ndim) ! resolution (CPPR), 0=max. character(len=12) :: type ! projectile type (shape) integer :: off(1:ndim) ! offset from center character(len=12) :: tprof ! Thermal profile flag character(len=12) :: pprof ! Porosity profile real*8 :: vel ! velocity real*8 :: angle(2) ! incidence angles (X/Z-dir, Y/Z-dir) real*8 :: vol ! volume real*8 :: rad(1:ndim) ! radius real*8 :: mid(1:ndim) ! coordinates of center real*8 :: sie ! specific internal energy real*8 :: temp ! temperature real*8 :: density ! density real*8 :: pressure ! pressure real*8 :: par_frac ! fraction of space occupied by particles real*8 :: par_range ! Particle size range factor about mean real*8 :: par_margin_h ! Particle margin on sides relative to particle size integer :: par_host_obj ! host object for particles character(len=12) :: par_dist ! Type of particle distribution end type type_obj type(type_obj), target, allocatable :: obj(:)
As you see, this datatype contains every information specified for each "object".
Please note that although in the input-file we distinguish between
iSALEalso layers are treated as
objectswith an infinite extent.
As you see, there are already some variables 'sie' (for specific internal energy), 'temp' (for temperature), 'density', and pressure which can be used to define the initial condition of each object. The only thing we need to do here is to add a new variable (let's call it "damage") into this structure:
type type_obj ... real*8 :: sie ! specific internal energy real*8 :: temp ! temperature real*8 :: density ! density real*8 :: pressure ! pressure real*8 :: damage ! damage ... end type type_obj
We save and close this file now.
Enable iSALE to read in the new variable¶
Now it's time to look how we can force
iSALE to read in the appropriate values from
asteroid.inp and store them into the 'damage' variable of the
obj-type. The library
ptools which is provided with
iSALE and automatically compiled and installed with
iSALE provides some helpful functions which ease the reading of input-parameters. We will show you now how to use them.
First, let's look into
Please look into the section '
! ------ Projectile parameters ---------------------------'. In this section all projectile-related parameters are read. Please look into the loop over all objects. There you find the following lines:
do i=1,obj_num oo=>obj(i) ... ! Optional: set projectile temperature, if not set, T_surf is used call ptool_get_dble("OBJTEMP",oo%temp,icol=i,defval=T_surf) ... enddo
The function '
ptool_get_dble' is provided with the
ptool-library. This line ensures that a double-precision-value (8-byte floating-point number) is read in from
asteroid.inp. The value is to be found in the line abbreviated with '
OBJTEMP' and in column 'icol=i' (i.e. the value for the first object is the first value, the value for the second object is the second value and so on). The value is read in and stored in
oo%temp (which is referring to
Please note that often pointers are defined to increase the readability of the code. The line
oo=>obj(i)means that the pointer
obj(i). This way you can use
oo%tempas an abbreviation of
The entry '
defval=T_surf' enables that
oo%temp is initialized with the surface-temperature (
T_surf) if this line of the value is not found in
According to this scheme, we now include some lines which read and initialize our damage-information for the projectile(s):
oo=>obj(i) ... ! Optional: set projectile temperature, if not set, T_surf is used call ptool_get_dble("OBJTEMP",oo%temp,icol=i,defval=T_surf) ! Optional: set initial damage-state of the projectile. ! If not set, define intact material (damage=0.) call ptool_get_dble("OBJDAM",oo%damage,icol=i,defval=0.D0) ...
With this line we ensure that the appropriate values are read in from the line beginning with "OBJDAM" and stored in obj(i)%damage. Per default, damage=0.D0 (intact material) is used, if the information is not specified in
Now, we also need to repeat these changes for the layer-information. Let's add the corresponding line(s) into the appropriate section:
do i=1,lay_num oo=>obj(obj_num+i) ... ! Define object properties as appropriate for layers oo%density = 0.D0; oo%sie = 0.D0; oo%pressure = P_surf; oo%temp = T_surf ! Optional: set initial damage-state of the projectile. ! If not set, define intact material (damage=0.) call ptool_get_dble("LAYDAM",oo%damage,icol=i,defval=0.D0) enddo
Once you add new lines into the source code, please ensure that you provide a good documentation for it. At least a small comment in the source code is quite helpful for other users.
Use new variable in the source code¶
Now it's time to use the stored values to influence the setup.
Let's save and close this file and open
This is part of the source-code of
iSALE-2D. In this routine, the cells are initialized with proper values prior to start the calculation. This is primarily done in the second loop (over all cells) which we will modify now. We add the following lines somewhere in the loop.
! initialize damage for non-void cells if (objID /= 0) damage(i,j) = obj(objID)%damage
Please note that these lines should be inserted after
objID is calculated and not within any other if-clause. A good place would be
do j = js,je !js,je do i = is,ie !is,ie ! Get cell's object number objID = cellType(i,j)%objID ... ... ! Get material ID of cell matID = cellType(i,j)%matID if (objID /= 0) then ! initialize damage for non-void cells damage(i,j) = obj(objID)%damage endif ! ----- ASSIGN MATERIAL PROPERTIES IN PROJECTILE(S) ----- ! if (objectFlag .and. matID /= VOID) then ... ... enddo enddo
With this line, we initialize the damage of the cell (
damage(i,j)) with the pre-defined initial damage value for the corresponding "object" (defined by
objID, which refers either to a projectile, a layer, or
For completeness, let us also output the object damage to the setup file along with all the other object properties. This is done in a subroutine called
output_model_setup_information at the end of the file. The relevant lines to change are:
... write(IOINIT,*) ' Object mass: ',oo%density*oo%vol write(IOINIT,*) ' Object porosity: ',mat(oo%mat)%alpha0 write(IOINIT,*) ' Object damage: ',oo%damage ...
Add parameter information to the database¶Now we are -- in principle -- ready to run the code. But before doing this, we should update the parameter database. This way, we give
iSALEPar, respectively) some meta-information about this parameter. This information is then used to
- automatically update the documentation (the
- prove automatically whether the values provided for this parameter are reasonable
- allow quick access to the meta-information from command-line (e.g.
iSALEPar --info LAYDAM)
- allow automatic cleaning-up or re-arranging of the input-file (e.g. grouping parameters by their meanings)
The database entries are found in
For easier handling of the meta-information, every parameter is described in it's own file, ending with
*.par. Let's use
GRIDSPC.paras a template and copy this file:
cp GRIDSPC.par LAYDAM.par
Now let's edit our newly created parameter file. The original file looks like this
<PARAM> ABBREV : GRIDSPC : DBLE DESC : grid spacing <INFO> The real-life spacing of cells in the high-resolution zone in metres. Two (or three) spacings can be defined: the first is the cell size in the horizontal direction; the second is the cell size in the vertical direction; the third is the cell size in the third dimension. If only one value is defined, the same spacing will be used for all coordinate directions. Note that when regridding (iSALE2D only), GRIDSPC is used to define the dimensions of the cell in the new grid; however, the only permissable values for GRIDSPC in this case are two times the GRIDSPC in the run </INFO> SECTION : MESH <CODE> DIM : 1 OPTIONAL : NO RANGE : 1.D-9 : RRANGE: 1.D-6 : 1.D6 COL: 1:3 </CODE> </PARAM>
Now we change some lines:
ABBREV : LAYDAM : DBLE
This way we define that the appropriate value is stored in the line beginning with "LAYDAM" of
DESC : initial layer damage
This line defines a short meaning of the parameter and is used e.g. in the manual or as a short description in
Next, we include some more comprehensive text explaining the meaning of this parameter. This text is found between
<INFO> The initial damage of each layer. By default, an initially intact material is assumed. Use this parameter to start with an arbitrary damage-value (between 0. and 1.). </INFO>
We now change the section from
TARGET. This way we assign this parameter to the group of "target-parameters". We also set
OPTIONAL : YES to define that this parameter is an optional one. To define the default-return value we include
DEFVAL : 0.D0.
We now adjust the physical (i.e. allowed) range for this parameter (
RANGE : 0.D0 : 1.D0) to assure that no values smaller than 0. and no values larger than 1.D0 are allowed to defined in the input-file for this parameter. The information for the recommended range (
RRANGE) can be removed, since this has the same values as the physical range.
The last entry (
COL) defines in which column(s) the value is found. Here we use
COL : 1 : (please don't forget the last ':') to state that an arbitrary number of columns might exist for this parameter (i.e. arbitrary number of layers).
Finally, the entire file will look as the following:
<PARAM> ABBREV : LAYDAM : DBLE DESC : initial layer damage <INFO> The initial damage of each layer. By default, an initially intact material is assumed. Use this parameter to start with an arbitrary damage-value (between 0. and 1.). </INFO> SECTION : TARGET <CODE> DIM : 1 OPTIONAL : YES DEFVAL : 0.D0 RANGE : 0.D0 : 1.D0 COL: 1 : </CODE> </PARAM>
In the same way we proceed now and insert the information for the projectile damage (
<PARAM> ABBREV : OBJDAM : DBLE DESC : initial projectile damage <INFO> The initial damage of each projectile. By default, an initially intact material is assumed. Use this parameter to start with an arbitrary damage-value (between 0. and 1.). </INFO> SECTION : OBJECT <CODE> DIM : 1 OPTIONAL : YES DEFVAL : 0.D0 RANGE : 0.D0 : 1.D0 COL: 1 : </CODE> </PARAM>
Build the code; verify the implementation and documentation¶
Now we build the code. To assure that all changes are really taken into account when compiling again, it is highly recommended to call
make clean before building again the code. This is especially true if you performed changes in a module, as we did here:
cd iSALE/trunk make clean make make install
Now we can run our first simulation to check whether everything works correctly. Therefore we change into your appropriate example-directory and call
Please note: You might also need to copy the new versions of the binary into this folder prior to starting the calculations.
If we compare our setup now with the original one (Chicxulub-example) we will see the differences if we visualize the damage-field:
Comparison of the original Chicxulub-setup (left) and the modified one (right). In each figure, the left half shows the material plot and the right half shows damage. As you see, our new setup matches our wishes: A totally damaged (damage=1; red color) sedimentary layer, a half damaged (damage=0.5; yellow color) granite layer and an intact (damage=0.; blue) mantle material.
Add simple test to check whether the input variable is read and used correctly¶
As any modification to a complex computer code may have unintended negative consequences, an important part of model development in large software projects is regression testing. For this reason, iSALE includes a suite of "tests" that should be run manually by developers prior to committing changes and that are run automatically by the iSALE buildbot whenever a change is committed to the repository.
It is best practice to create a new test for every modification made to iSALE. A tutorial describing how to create a test is provided here. However, as time is limited in this tutorial, we will modify an existing test to ensure that the new feature (prescribing an initial object/layer damage) works, rather than creating a new test.
One of iSALE's tests,
options_check_2D, verifies that for each example problem the input files create no errors and that the simulations are set-up correctly. We can run this test individually as follows (from the iSALE root-directory):
cd iSALE/tests/options_check_2D ../testharness.py -f options_check_2D.xml
testharness.py is the program that performs the test; the xml file prescribes the definition of the test.
The test should pass at this stage, because we have not made any changes to it and no example uses the new feature. To include the new feature in this test, let us modify the Chicxulub example so that the impactor is assumed to be fully damaged. First open the relevant example input file (in its repository location):
Then, add the new object property (and save the file):
------------------- Projectile ("Object") Parameters -------------------- OBJNUM number of objects : 1 OBJRESH CPPR horizontal : 18 OBJVEL object velocity : -1.2D4 OBJMAT object material : granite OBJTYPE object type : SPHEROID OBJDAM object damage : 1.D0
Now we can re-run the test:
../testharness.py -f options_check_2D.xml
If the new feature is successfully implemented, the test should still pass. However, the test is not very useful as it does not specifically verify that the impactor has indeed been assigned a damage of 1. A simple way to do this is to add a simple additional pass_test to the xml file:
The bits relevant to the Chicxulub example are, the variable definitions:
<!--Chicxulub example--> <variable name="ChicxSetupFile" language="python"> from isale_tools import read_contents_of_file ChicxSetupFile = read_contents_of_file("./Chicxulub/INFO/setup_report.txt") </variable> <variable name="ChicxErrorFile" language="python"> from isale_tools import read_contents_of_file ChicxErrorFile=read_contents_of_file("./Chicxulub/errors.txt0","Error file does not exist. Probably because of invalid input parameter(s).") </variable> <variable name="ChicxOutputFile" language="python"> from isale_tools import read_contents_of_file ChicxOutputFile=read_contents_of_file("./Chicxulub/output.txt0") </variable>
Which create three "variables" containing the contents of three files: the setup report file, the error file and the output file.
And the tests:
<!--Chicxulub example--> <test name="nothing in error file" language="python"> import re print "Chicxulub test" print ChicxulubErrorFile assert(not(re.search(" ",ChicxErrorFile,re.M))) </test> <test name="model finished set-up correctly" language="python"> import re assert(re.search("SETTINGS FINISHED, STOPPING.",ChicxSetupFile,re.M)) </test> <test name="model terminated correctly" language="python"> import re assert(re.search("END OF SIMULATION",ChicxOutputFile,re.M)) </test>
The first test passes if no text is found in the error file. The second test passes if the string
SETTINGS FINISHED, STOPPING. is found in the setup file. The third test passes if the string:
END OF SIMULATION is found in the output file.
As the setup file now includes a line defining the damage of the object, we can simply add one additional test that verifies that the object damage is 1, by searching the setup file for the relevant string. I.e., we add the following test to the three above...
<test name="check whether object damage is correct" language="python"> import re assert(re.search("Object damage: 1.00",ChicxSetupFile,re.M)) </test>
Run the tests to verify that changes are good¶
Now that we have a useful test, we should run this test and all the other tests to verify that the changes we have made have not had any unintended consequences. Assuming that we have gone back to the iSALE root directory:
Commit changes to repository¶
Assuming that all of the tests pass, we are ready to commit our changes to the repository. Before we do this, we need to make sure that our version of the code is up-to-date as we cannot commit changes to an old revision. So first we update:
All being well, we will already be at the most recent revision, but if some has modified the code since our last update watch out for any conflicts and take the necessary steps to resolve these. For more information about resolving conflicts, see: this page.
Next, we should check to see which files have been modified. This can easily be accomplished with:
svn status -q
Which informs us that the following files have changed:
M iSALE/tests/options_check_2D/options_check_2D.xml M iSALE/src/2D/setup_cells.F90 M iSALE/src/shared/mod_setup_params.F90 M iSALE/src/shared/read_asteroid.F90 M iSALE/example/Chicxulub/asteroid.inp
Notice that the parameter files that we have created:
OBJDAM.par, are not listed. This is because the
-q flag omits files not already in the repository and we have not yet added these files to the repository---it does not know that they exist. To add these files we:
svn add iSALE/parameters/OBJDAM.par iSALE/parameters/LAYDAM.par
Now if we
svn status -q again, we should see listed all the files that we have changed or added:
M iSALE/tests/options_check_2D/options_check_2D.xml A iSALE/parameters/LAYDAM.par A iSALE/parameters/OBJDAM.par M iSALE/src/2D/setup_cells.F90 M iSALE/src/shared/mod_setup_params.F90 M iSALE/src/shared/read_asteroid.F90 M iSALE/example/Chicxulub/asteroid.inp
At this stage we are almost ready to commit our modifications; we just need to compose a suitably informative message for other developers to accompany the commit. For large commits it is best to compose this in a text editor and save the message in a temporary text file. But for our purposes a short note can be added with the
svn ci -m"Added new object and layer parameter to define the initial damage. This parameter is now used in iSALE2D but has not yet been used in iSALE3D. Parameter descriptions for the two new parameters have been added; the Chicxulub example has been modified to use a fully damaged impactor and the options_check_2D test has been improved to verify that object damage is set correctly in this example."
Note that we have not specified which files to commit, so all modified files will be committed, including those in subdirectories. Commits can be made more specific by listing each file individually.