make_singly_periodic_grid.sh

#!/usr/bin/env bash
# Phillip Wolfram 02/04/2015 from https://gist.github.com/pwolfram/8724b45d788bcdfa20f7
# stored in gist: 76e9bd312c7a3e048f32: 
function update_gist { #{{{
  echo 'updating the gist'
  gist  make_singly_periodic_grid.sh -u 76e9bd312c7a3e048f32
  echo 'done'
} #}}}
update_gist

############################################################
# set up test case definitions
############################################################
function use_250m_res { #{{{

  echo 'building high resolution grid'
  export Nhoriz=100
}
function use_highest_res { #{{{

  echo 'building high resolution grid'
  export Nhoriz=400
}
#}}}
function use_high_res { #{{{

  echo 'building high resolution grid'
  export Nhoriz=200
}
#}}}
function use_med_res { #{{{

  echo 'building medium resolution grid'
  export Nhoriz=100
}
#}}}
function use_low_res { #{{{

  echo 'building low resolution grid'
  export Nhoriz=50
}
function use_ultra_low_res { #{{{

  echo 'building lowest resolution grid'
  export Nhoriz=10
}
#}}}
#export W=1000000.000
#export L=2000000.000
export W=1000.000
export L=2000.000
export W=25000
export L=50000
export bottom=2500.0
export HOMEDIR=`pwd`
export METIS='gpmetis'
export NPROCS=6
echo 'currently in '$HOMEDIR' with '$NPROCS' processors'

############################################################
# initialize basic grid parameters
############################################################
#use_low_res
#use_med_res
#use_250m_res 
use_ultra_low_res
export dc=`echo 'print '$W'/'$Nhoriz | xargs -I {} python -c {}`
export Nvert=`echo 'from math import sqrt; import numpy as np; n = int(np.round(2.0/sqrt(3.0)*'$L'/'$dc')); n -= np.mod(n,2); print n' | xargs -I {} python -c {}`

function download_code { #{{{

  # periodic hex generator
  git clone git@github.com:pwolfram/MPAS-Tools.git periodic_hex_generator
 
  # mesh converter tools
  git clone git@github.com:pwolfram/MPAS-Tools.git mesh_converter_culler
  cd $HOMEDIR

  # init core (new to setup test case)
  git clone git@github.com:pwolfram/MPAS.git MPAS_tracerInfrastructure
  cd MPAS_tracerInfrastructure
  git co ocean/ZISOcase
  cd $HOMEDIR

}
#}}}
function compile_code { #{{{

  # mesh converter / cell culer
  make -C mesh_converter_culler/grid_gen/mesh_conversion_tools
  # periodic hex generator
  make -C  periodic_hex_generator/grid_gen/periodic_hex/
  # init core (new)
  make -C MPAS_tracerInfrastructure gfortran CORE=ocean AUTOCLEAN=true DEBUG=false
  # mpas draw
  #sed -i .back \
  #  -e 's/PLATFORM=_LINUX.*/#PLATFORM=_LINUX/g' \
  #  periodic_hex_generator/visualization/mpas_draw/Makefile
  #make -C periodic_hex_generator/visualization/mpas_draw/ 

}
#}}}
function clean_up_code { #{{{
  rmtrash MPAS_init_core/ 
  rmtrash MPAS_tracerInfrastructure/ 
  rmtrash mesh_converter_culler/ 
  rmtrash periodic_hex_generator/ 
}
#}}}
function clean_up_intermediates { #{{{
  rmtrash *.nc 
  rmtrash step*nc 
  rmtrash *log* 
  rmtrash *info* 
  rmtrash *graph* 
  rmtrash namelist.ocean.init.*  
}
#}}}
function test_dimensions { #{{{
  python test_ranges.py -f $1
}
#}}}
############################################################
# step 0: preliminaries 
############################################################
clean_up_code
clean_up_intermediates
download_code
compile_code

############################################################
# step 1: creating an arbitrary grid
############################################################
function build_grid { #{{{

  # script to build up a singly-periodic MPAS-O mesh
  cd periodic_hex_generator/grid_gen/periodic_hex/

  # update the namelist
  sed -i .back \
    -e 's/dc =.*/dc = '$dc',/' \
    -e 's/nx =.*/nx = '$Nhoriz',/' \
    -e 's/ny =.*/ny = '$Nvert',/' \
  namelist.input

  ./periodic_grid &> step1.log


  cp grid.nc $HOMEDIR/step1.nc

  cd $HOMEDIR

  #makes step1.nc
} #}}}
build_grid

############################################################
# step 2: Converting arbitrary grid to an MPAS mesh
############################################################
function make_mpas_mesh { #{{{
  
  # use mesh converter
  ./mesh_converter_culler/grid_gen/mesh_conversion_tools/MpasMeshConverter.x step1.nc step2.nc &> step2a.log
 
} #}}}
make_mpas_mesh

#############################################################
## step 3: identifying cells for removal
#############################################################
function make_singly_periodic { #{{{

  ## use init core to make mesh singly periodic
  #cp MPAS_init_core/namelist.ocean_init.defaults namelist.ocean_init
  #sed -i .back \
  #  -e 's/config_do_restart =.*/config_do_restart = .false./' \
  #  -e 's/config_vert_levels =.*/config_vert_levels = '1'/' \
  #  -e 's/config_input_name =.*/config_input_name = "step2.nc"/' \
  #  -e 's/config_output_name =.*/config_output_name = "step3.nc"/' \
  #  -e 's/config_write_cull_cell_mask =.*/config_write_cull_cell_mask = .true./' \
  #  -e 's/config_remove_y_boundary =.*/config_remove_y_boundary = .true./' \
  #  namelist.ocean_init

  #$METIS graph.info $NPROCS &> step3_metis.log
  #mpirun -np $NPROCS ./MPAS_init_core/ocean_model &> step3_mpas.log

  ## take care of periodicity
  #python add_periodic_flags.py step3.0000-01-01_00.00.00.nc $L `echo $L*10 | bc -l`

  # mark y-periodic cells to cull
  cp step2.nc step3.nc
  #./periodic_hex_generator/grid_gen/periodic_hex/mark_periodic_boundaries_for_culling.py -f step3.nc --yonly
  ./periodic_hex_generator/grid_gen/periodic_hex/mark_periodic_boundaries_for_culling.py -f step3.nc 

} #}}}
make_singly_periodic

############################################################
# step 4: removal of cells
############################################################
function cell_removal { #{{{

  # note: could be mitigated by test case itself

  # remove culled cells
  ./mesh_converter_culler/grid_gen/mesh_conversion_tools/MpasCellCuller.x step3.nc step4.nc &> step4a.log
  
  ./mesh_converter_culler/grid_gen/mesh_conversion_tools/MpasMeshConverter.x step4.nc &> step4b.log
  mv mesh.nc step4.nc
  #mv culled_graph.info graph.info
  
  ## take care of periodicity
  #python add_periodic_flags.py step4.nc $W `echo $L*10 | bc -l`

  # visualize
  #~/Documents/MPAS-Tools/visualization/mpas_draw/MpasDraw.x step4.nc

} #}}}
cell_removal
test_dimensions step4.nc

############################################################
# step 5: applying initial conditions
############################################################
function apply_initial_conditions { #{{{
  
  #sed -i .back2 \
  #  -e 's/config_do_restart =.*/config_do_restart = .false./' \
  #  -e 's/config_vert_levels =.*/config_vert_levels = '-1'/' \
  #  -e 's/config_input_name =.*/config_input_name = "step4.nc"/' \
  #  -e 's/config_output_name =.*/config_output_name = "step5.nc"/' \
  #  -e 's/config_write_cull_cell_mask =.*/config_write_cull_cell_mask = .false./' \
  #  -e 's/config_remove_y_boundary =.*/config_remove_y_boundary = .false./' \
  #  -e 's/config_initial_conditions =.*/config_initial_conditions = "ziso"/' \
  #  -e 's/config_iso_zonal_periodic_north_wall_y =.*/config_iso_zonal_periodic_north_wall_y =   '$W'./' \
  #  -e 's/config_iso_zonal_periodic_main_channel_depth =.*/config_iso_zonal_periodic_main_channel_depth = '$bottom'/' \
  #  namelist.ocean_init

  cp MPAS_tracerInfrastructure/defaults/namelist.ocean.init.ziso namelist.ocean
  # namelist update
  sed -i .back2 \
    -e "s/config_ocean_run_mode = .*/config_ocean_run_mode = 'init'/" \
    -e "s/config_init_configuration = .*/config_init_configuration = 'ziso'/" \
    -e "s/config_vert_levels =.*/config_vert_levels = 100/" \
    -e "s/config_write_cull_cell_mask =.*/config_write_cull_cell_mask = false/" \
    -e "s/config_ziso_use_slopping_bathymetry =.*/config_ziso_use_slopping_bathymetry = false/" \
    -e "s/config_use_bulk_wind_stress =.*/config_use_bulk_wind_stress = true/" \
    -e "s/config_ziso_use_slopping_bathymetry =.*/config_ziso_use_slopping_bathymetry = false/" \
    -e "s/config_use_activeTracers_surface_restoring =.*/config_use_activeTracers_surface_restoring = true/" \
    -e "s/config_use_activeTracers_interior_restoring =.*/config_use_activeTracers_interior_restoring = true/" \
    namelist.ocean

  # streams update
  python - <<END
from namelist_interface import XMLList
stream = XMLList('streams.ocean')
stream.modify('input_init', 'filename_template', 'step4.nc')
stream.modify('output_init', 'filename_template', 'step5.nc')
stream.write()
END

 
  # run init core
  $METIS graph.info $NPROCS &> step5_metis.log
  mpirun -np $NPROCS ./MPAS_tracerInfrastructure/ocean_model &> step5_mpas.log

  ## take care of periodicity
  #python add_periodic_flags.py step5.0000-01-01_00.00.00.nc $W `echo $L*10 | bc -l`

} #}}}
#apply_initial_conditions
#test_dimensions step5*.nc
#
############################################################
# step 6: map coordinates to the cylinder
############################################################
function map_to_cylinder { #{{{
  
  cp step5*nc step6.nc
  # modify in-place with a python script
  gist  map_to_cylinder.py -u 76e9bd312c7a3e048f32
  python map_to_cylinder.py -f step6.nc -w $W
} #}}}
#map_to_cylinder
############################################################
# step 7: run
############################################################
# mv step5.nc to init.nc and mesh.nc
# mv init_forcing_data.nc to forcing_data.nc
# edit namelist and stream files from defaults
# run the test case

map_to_cylinder.py

#!/usr/bin/env python
# Phillip J. Wolfram 02/05/2015 
# map points in grid onto a cylinder and replace values

import netCDF4
import numpy as np
import copy

def remap_points(f_in, name, width):
    # get the data 
    xp = copy.deepcopy(f_in.variables['x'+name][:])
    yp = copy.deepcopy(f_in.variables['y'+name][:])
    zp = copy.deepcopy(f_in.variables['z'+name][:])

    # compute the coordinate transform
    twopi = 2.0*np.pi
    radius = width/twopi
    theta = xp[:]/width
    xc = radius*np.sin(twopi*theta)
    yc = radius*np.cos(twopi*theta)
    zc = yp[:] - np.mean(yp[:])

    # place new transformed coordinates back in the file
    f_in.variables['x'+name][:] = xc[:]
    f_in.variables['y'+name][:] = yc[:]
    f_in.variables['z'+name][:] = zc[:]

def compute_lat_lon(f_in, name):
    pass

def remap_to_cylinder(f_in, width):
    f_in = netCDF4.Dataset(f_in, 'r+')

    #compute_lat_lon(f_in, 'Cell')
    #compute_lat_lon(f_in, 'Vertex')
    #compute_lat_lon(f_in, 'Edge')

    remap_points(f_in, 'Cell', width)
    remap_points(f_in, 'Vertex', width)
    remap_points(f_in, 'Edge', width)

    f_in.close()

if __name__ == "__main__":
    from optparse import OptionParser

    # Get command line parameters
    parser = OptionParser()
    parser.add_option("-f", "--file", dest="f_in",
                      help="file to open for testing", \
                      metavar="FILE")
    parser.add_option("-w", "--width_of_plane", dest="width",
                      help="width of plane's x coordinate in m", \
                      metavar="REAL")

    options, args = parser.parse_args()

    if not options.f_in:
        parser.error("Filename is a required input.")
    if not options.width:
        parser.error("Width is a required input.")

    remap_to_cylinder(options.f_in, float(options.width))