Skip navigation
NASA Logo, National Aeronautics and Space Administration

Thoughts on Modeling with TDD

2 Posts tagged with the pfunit tag
0

When an unsupported situation (e.g. illegal value) arises in a scientific model, the usual implementation is to stop the model.   In more sophisticated cases, some possibly-useful message about why the model stopped is generated.  However, the logic for handling these rare cases is often buggy - precisely because the situation in rarely encountered!   At any rate, this report-and-abort approach is problematic for TDD because a test which encounters the unsupported situation fails to complete, as opposed to returning with an error condition.   (Note: frameworks in languages such as python and java can exit gracefully from these types of situations, but Fortran has no true exception handing mechanism.)

 

But there is another path.

 

If instead of halting the execution, the routine in question were to "throw" an exception, the framework would be back-in-business.   Even though Fortran does not support true exceptions, the various Fortran unit testing frameworks (and pFUnit in particular) provide means to throw pseudo-exceptions that can then be handled at some other level in the software.   (For the curious, the bit of exception handling that is missing is the ability to pop back up multiple levels of the procedure stack.)

 

So how do we drive this approach with TDD?  Simple.  We write a test which intentionally creates the "unsupported" condition, and then checks that the desired exception is thrown.  "Unsupported" is in quotes, because technically we are now supporting it, albeit in a limited fashion.    In order to pass this test, the implementation must throw the desired exception, but only under the appropriate circumstances.  (Otherwise the other more-typical tests would then signal failure due to the unexpected exception.)

 

The frustrating downside of this approach is that the source code (as opposed to the test code) must now link to the testing framework.   One can use CPP, to minimize the impact on the other developers, but CPP is its own evil.

0

Nearly all of the modeling of physical phenomena is based on conservation principles with the corresponding conservation equations represented by partial (or ordinary) differential equations. These equations are almost always impossible to solve analytically and thus one must immediately start choosing  numerical techniques appropriate for each problem.

 

In our first post we take a bird's eye view of modeling and choose the decay problem as a first example to demonstrate the use of TDD. In mathematical terms, this problem can be expressed as as an initial value problem:

 

dy(t)/dt = -ky, y(0)=y0    (1)

 

where y is concentration of a substance, t is time, and k is a positive constant parameter. Note that with k=0 there will be no change and y remains unchanged at its initial value. We also note that y will gradually decrease with time at a rate in proportion to concentration itself at any time instance.

 

So, how do we proceed in creating our first first numerical prediction model using TDD? What should be our first unit test? What should we test?

 

First we must decide what numerical method method to use. Let's start with the Euler method. In that case (1) can be discretized as follows:

 

y(n+1) = y(n)(1 + ht)     (2)

 

where h is the time step and n is the time index. We know from (1) that for k=0  y remains unchanged at its initial value so (2) must satisfy that condition. Let's create a module to contain our test:

 

Note: In this and all future posts we will be mainly making use of the fortran programming language and will be performing our unit tests under the pFUnit testing framework. For more information please see Modeling with TDD

 

module testIntegrators_mod
   use pFunit                ! unit testing framework
   use Integrators_mod       ! our code implememtation
   implicit none
   private
   public testEuler
   integer, parameter :: dp = selected_real_kind(14) ! 64 bit precision

contains

   subroutine testEuler
   implicit none
   ! y is the solution to the decay equation    
   real(kind=dp), allocatable, dimension(:) :: y
   real(kind=dp) :: y0       ! initial condition
   real(kind=dp) :: timeStep ! in seconds
   real(kind=dp) :: kappa    ! in seconds^-1
   integer :: numSteps
   integer :: i

! test condition that if k=0, yf=y0

   timeStep = 3600
   kappa = 0.0
   numSteps = 86400/timeStep
   allocate(y(numSteps))
   y0 = 1.0
   call Euler(y, y0, timeStep, numSteps)
   call assertEqual(y0, y(numSteps), tolerance=1e-8)
   deallocate(y)

   end subroutine testEuler

end module testIntegrators_mod

 

Note that our code implementation, Integrators_mod, does not yet exist and as such our code test module will not even compile! This is what TDD intends: let the unit testing drive the code development. However even at this point in our simple exercise several questions come to mind: how simple (or complex) should our first test be? what assumptions should be make as we write the test, say regarding  timeStep and numSteps? To what tolerance should be expect our result to be true? There is probably no clear-cut answer.

 

So now we must write the simplest unit of code that will make our test pass. Again, what constitues the simplest test? Here are two possibilities. First,

 

Module Integrators_mod
   implicit none
   private
   public Euler
   integer, parameter :: dp = selected_real_kind(14) ! 64 bit precision

contains

   subroutine Euler(yf, y0, timeStep, numSteps, kappa)
   real(kind=dp), intent(out) :: yf(:)
   real(kind=dp), intent(in)  :: timeStep
   real(kind=dp), intent(in)  :: y0
   real(kind=dp), intent(in)  :: kappa
   integer, intent(in)        :: numSteps

   integer :: i
   real :: time
   real :: yy0

   yy0 = y0
   do i = 1,numSteps
      time = time + timeStep
      yf(i) = yy0*(1.0 - timeStep*kappa)
      yy0 = yf(i)
   end do

   end subroutine Euler

end Module Integrators_mod

 

Second, we note that subroutine Euler could be as simple as

 

   subroutine Euler(yf, y0, timeStep, numSteps, kappa)
   real(kind=dp), intent(out) :: yf(:)
   real(kind=dp), intent(in)  :: timeStep
   real(kind=dp), intent(in)  :: y0
   real(kind=dp), intent(in)  :: kappa
   integer, intent(in)        :: numSteps
   yf(numSteps) = y0
   end subroutine Euler
 

 

and this will pass the test. What is the best approach? To be continued.



USAGov logo NASA Logo - nasa.gov