Team 3 Solution

From Interactive System for Ice sheet Simulation
Jump to: navigation, search
  • FORTRAN mountain glacier

Result for bed slope 0.025 Kees result.png

Compile the program with gfortran -o simple_ice_model simple_ice_model_modules.f90 simple_ice_model.f90. The python script does the compiling automatically, so just running python will do it all.

Source: Kees_g3.tar.gz‎


program simple_ice_model
  use simple_ice_model_modules
  implicit none
  ! local variables
  real(8), parameter    :: pi = 3.1415926536
  integer, parameter :: grid=51                    ! Number of nodes
  real(8),    parameter :: dt = 0.1                ! length time step 
  integer, parameter :: t_final=1000                ! integration stop time
  integer :: t_step, nt                                ! current time step, number of time steps
  real(8), parameter :: xl=0                       ! start of domain
  real(8), parameter :: xr=50e3                    ! end of domain 
  real(8) :: dx                                    ! node spacing
  real(8) :: time                                  ! current time
  real(8) :: const                                 ! a constant
  !! unstaggered grid
  real(8), dimension(:), allocatable :: xx         ! node positions
  real(8), dimension(:), allocatable :: thick      ! ice thickness
  real(8), dimension(:), allocatable :: dthick_dt  ! ice thickness time derivative
  real(8), dimension(:), allocatable :: mass_b     ! mass balance
  real(8), dimension(:), allocatable :: bed        ! bed elevation
  !! staggered grid
  real(8), dimension(:), allocatable :: diffus     ! diffusion
  real(8), dimension(:), allocatable :: xx_s       ! xx on staggered
  !! physical constants
  real(8), parameter :: rho=920, g_grav=9.8, A_rate=1e-16, n_glen=3.
  real(8), parameter :: dbed_dx=0.025, mass_b_0=4., mass_b_1=0.2e-3 !xxx
  !! BC
  real(8), parameter :: thickl=0, thickr=0 !xxx
  !! counters and such
  integer :: ii, jj
  integer :: errstat                            ! for error checking
  !! functions
  ! let us allocate some memory
  call checkerr(errstat,"failed to allocate xx")
  call checkerr(errstat,"failed to allocate thick")
  call checkerr(errstat,"failed to allocate dthick_dt")
  call checkerr(errstat,"failed to allocate mass_b")
  call checkerr(errstat,"failed to allocate bed")
  call checkerr(errstat,"failed to allocate diffus")
  call checkerr(errstat,"failed to allocate xx_s")
  !! setup
  dx = (xr-xl)/(grid-1)
  do ii=1,grid
     xx(ii) = xl + dx*(ii-1)
     thick(ii) = 0   ! IC
     mass_b(ii) = mass_b_0 - mass_b_1*xx(ii)
     bed(ii) = -dbed_dx * xx(ii)
  end do
  const = 2*A_rate/(n_glen+2)*(rho*g_grav)**n_glen
  nt = floor(t_final/dt) ! number of time steps
  !! output, use a pipe to get into a file
  print * , 'x'
  print * , xx
  print * , 'mass_b'
  print * , mass_b
  print * , 'bed'
  print * , bed
  print *, 'thick'
  do ii=1,nt
     ! compute the diffusivity: diffus
     diffus = calc_diffus(thick, bed, dx, const, n_glen)
     ! assemble dthick_dt
     dthick_dt = pad(diff_1(diffus*diff_1((thick+bed), dx), dx), thickl, thickr) + mass_b
     ! time step
     thick = time_step(thick, dthick_dt, dt)
     ! BC
     thick(1) = thickl
     do jj=1,grid
        if (thick(jj)<0.) then
          thick(jj) = 0.
        end if
     end do
     print *, thick
  end do
  subroutine checkerr(errstat,msg)
    !! the error checker for allocation
    implicit none
    integer,      intent(in) :: errstat
    character(*), intent(in) :: msg 
    if (errstat /= 0) then
       write(*,*) "ERROR:", msg
    end if
  end subroutine checkerr
end program simple_ice_model


module simple_ice_model_modules
  function on_stagger(xx)
    !! interpolates a vector on the normal grid to the staggered grid
    implicit none
    real(8), intent(in), dimension(:):: xx
    real(8), dimension(size(xx)-1):: on_stagger
    integer :: ii
    do ii=1,(size(xx)-1)
       on_stagger(ii) = (xx(ii+1)+xx(ii))/2
    end do
  end function on_stagger
  function calc_diffus(thick, bed, dx, const, n_glen)
    ! calculates the nonlinear diffusivity
    implicit none
    real(8), intent(in), dimension(:):: thick, bed
    real(8), dimension(size(thick)-1):: calc_diffus
    real(8), intent(in):: const, n_glen, dx
    integer :: ii
    calc_diffus = const * (on_stagger(thick))**(n_glen+2) * (diff_1((thick+bed), dx))**(n_glen-1)
  end function calc_diffus
  function diff_1(yy, dx)
    ! calculates the first staggered derivative of a normal vector 
    implicit none
    real(8), intent(in), dimension(:):: yy
    real(8), intent(in) :: dx
    real(8), dimension(size(yy)-1):: diff_1
    integer :: ii
    do ii=1,(size(yy)-1)
       diff_1(ii) = (yy(ii+1)-yy(ii))/dx
    end do
  end function diff_1
  function time_step(yy, dyy_dt, dt)
    ! does a euler time step
    implicit none
    real(8), intent(in), dimension(:):: yy, dyy_dt
    real(8), intent(in) :: dt
    real(8), dimension(size(yy)):: time_step
    time_step = yy + dt*dyy_dt
  end function time_step
  function pad(on_stag, left, right)
    ! pad a vector with left and right, used after double staggering
    implicit none
    real(8), intent(in), dimension(:):: on_stag
    real(8), intent(in) :: left, right
    real(8), dimension(size(on_stag)+2):: pad
    integer :: ii
    pad(1) = left
    pad(size(on_stag)+2) = right
    do ii=1,size(on_stag)
       pad(ii+1) = on_stag(ii)
    end do
  end function pad
end module simple_ice_model_modules

Python runner and plotter:

The runner and plotter. To make it work the fortran files need to be named as above

import numpy as N
import pylab as P
import os
from py_fortran_tools import *
data_output_file = 'tmp.dat'
f_output_file = 'simple_ice_model'
f_input_files = ['simple_ice_model_modules.f90', 'simple_ice_model.f90']
compile_and_run(f_output_file, f_input_files, data_output_file)
# get the data
dict_ = parse_output(data_output_file)
# plot it
shape_ = dict_['thick'].shape
P.plot(dict_['x']/1e3, dict_['thick'][-1,:]+dict_['bed'], label='At $t_{final}$')
P.plot(dict_['x']/1e3, dict_['thick'][N.floor(shape_[0]/2),:]+dict_['bed']
       , label='at $t_{final}$/2')
P.plot(dict_['x']/1e3, dict_['thick'][-100,:]+dict_['bed'],
       label='100 time steps before $t_{final}$')
P.xlabel('x (km)')
P.ylabel('b+H (m)')

Using the module

from __future__ import with_statement
import numpy as N
import pylab as P
import os
def parse_output(data_output_file):
    This function parses the output of a fortran program generated by
    print *, 'x'
    print *, x
    input: data_output_file
    output: dictionary with keys corresponding to the variable
    # parse the file
    with open(data_output_file, 'r') as fil:
        dict_ = {}
        # go through all the lines in the file
        for line in fil:
            ll = line.strip()
            if ll[0].isalpha():
                # if the first character in a line is a letter
                # it's the beginning of a new variable
                key = ll
                # make a new dict_ entry and initialise as empty list
                dict_[key] = []
                # otherwise read the data
                ll = ll.split() # split at white space
                num = [float(el) for el in ll] # convert to float
                dict_[key].append(num) # append to the list
        # convert to numpy arrays
        for key in dict_:
            dict_[key] = N.array(dict_[key]).squeeze()
    return dict_
def compile_and_run(f_output_file, f_input_files, data_output_file):
    Compiles a fortran program using
     - f_output_file as the executable produced
     - f_input_files as a list of the source code files (modules first)
     - data_output_file as the data file into which the output to
                        standard output of the fortran program
                        is piped
    # compile the fortran program
    command1 = 'gfortran -o %s %s' %(f_output_file, ' '.join(f_input_files))
    exit_ = os.system(command1)
    # some primitve error checking
    if exit_!=0:
        raise(TypeError('\nCompliling with command "%s" failed' %command1))
        print 'Compiling ... ok'
    # run it
    command2 = './%s > %s' %(f_output_file,  data_output_file)
    exit_ = os.system(command2)
    # some primitve error checking
    if exit_!=0:
        raise(TypeError('\nRunning of fortran program with command "%s" failed.' %command2))
        print 'Running ... finished'