Project

General

Profile

Tutorial: How to add a new parameter

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.

Currently, 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 asteroid.inp

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 install will 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 asteroid.inp:

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 asteroid.inp:

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

iSALE/trunk/iSALE/src/shared/mod_setup_params.F90

In this file, you find a so-called "derived datatype" type_obj:

  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 layers and objects (=projectiles), within iSALE also layers are treated as objects with 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

iSALE/trunk/iSALE/src/shared/read_asteroid.F90

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 obj(i)%temp).

Please note that often pointers are defined to increase the readability of the code. The line oo=>obj(i) means that the pointer oo refers to obj(i). This way you can use oo%temp as an abbreviation of obj(i)%temp.

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 asteroid.inp.

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 asteroid.inp.

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

iSALE/trunk/iSALE/src/2D/setup_cells.F90

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 VOID, if objID==0).

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 iSALE (or iSALEPar, respectively) some meta-information about this parameter. This information is then used to
  1. automatically update the documentation (the iSALE-manual)
  2. prove automatically whether the values provided for this parameter are reasonable
  3. allow quick access to the meta-information from command-line (e.g. iSALEPar --info LAYDAM)
  4. allow automatic cleaning-up or re-arranging of the input-file (e.g. grouping parameters by their meanings)

The database entries are found in

cd iSALE/trunk/iSALE/parameters

For easier handling of the meta-information, every parameter is described in it's own file, ending with *.par. Let's use GRIDSPC.par as 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 asteroid.inp.

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 asteroid.inp (2nd column).

Next, we include some more comprehensive text explaining the meaning of this parameter. This text is found between <INFO> and </INFO>:

<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 MESH to 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 (OBJDAM.par):

<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

iSALE2D -C

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:

Chicxulub scenario: Comparison between original (left) and modified (right) setup.

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):

nano ../../example/Chicxulub/asteroid.inp

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:

nano options_check_2D.xml

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:

make tests

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:

svn 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: LAYDAM.par and 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 -m flag:

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.

Chicxulub_compare_original_damage.png View - Chicxulub scenario: Comparison between original (left) and modified (right) setup. (4.43 KB) Dirk Elbeshausen, 07/25/2013 04:09 PM