#!/usr/bin/env python
#########################################################
#   
#  Main file for parallel mesh testing.
# 
#  This is a modification of the run_parallel_advection.py
# file.
#
#  *) The test files currently avaliable are of the form
# test*.out, eg test_5l_4c.out. The term infront of the l
# corresponds to the number of levels of refinement
# required to build the grid, i.e. a higher number
# corresponds to a finer grid. The term infront of the c
# corresponds to the number of processors. 
#
# *) The (new) files that have been added to manage the
# grid partitioning are
#    +) mg2ga.py: read in the test files.
#    +) build_submesh.py: build the submeshes on the host
# processor.
#    +) build_local.py: build the GA mesh datastructure
# on each processor.
#    +) build_commun.py: handle the communication between
# the host and processors
#
# *) Things still to do:
#    +) Overlap the communication and computation: The
# communication routines in build_commun.py should be
# interdispersed in the build_submesh.py and build_local.py
# files. This will overlap the communication and
# computation and will be far more efficient. This should
# be done after more testing and there more confidence in
# the subpartioning.
#    +) Much more testing especially with large numbers of
# processors.
#  Authors: Linda Stals, Steve Roberts and Matthew Hardy,
# June 2005
#
#
# 
#########################################################

import pypar    # The Python-MPI interface
import sys
from os import sep
sys.path.append('..'+sep+'pyvolution')

# advection routines

from config import g, epsilon
from Numeric import allclose, array, zeros, ones, Float
from parallel_advection import *
from Numeric import array

# mesh partition routines

from mg2ga import *
from build_submesh import *
from build_local import *
from build_commun import *
import pdb

# define the initial time step

class Set_Stage:
    """Set an initial condition with constant water height, for x<x0
    """

    def __init__(self, x0=0.25, x1=0.5, h=1.0):
        self.x0 = x0
        self.x1 = x1
        self.h  = h

    def __call__(self, x, y):
        return self.h*((x>self.x0)&(x<self.x1))

# read in the processor information

numprocs = pypar.size()
myid = pypar.rank()
processor_name = pypar.Get_processor_name()

print "PROCESSOR ", myid, " HAS STARTED"

# if this is the host processor

print 'trace'
pdb.set_trace()
print 'after trace'


if myid == 0:

    # read in the test files
    
    f=open('test_3l_2c.out', 'r')
    [nodes, triangles, boundary, triangles_per_proc] = mg2ga(f)
    quantities = {}
    
    # subdivide the mesh
    
    submesh = build_submesh(nodes, triangles, boundary, quantities, \
                            triangles_per_proc)

    # send the mesh partition to the appropriate processor
    
    for p in range(1, numprocs):
      send_submesh(submesh, triangles_per_proc, p)

    hostmesh = extract_hostmesh(submesh)
    [points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict] = \
             build_local_mesh(hostmesh, 0, triangles_per_proc[0], \
                              numprocs)
    
# read in the mesh partition that belongs to this
# processor (note that the information is in the
# correct form for the GA data structure

else:
    [points, vertices, boundary, quantities, ghost_recv_dict, full_send_dict] = rec_submesh(0)

# define the computation domain

domain = Parallel_Domain(points, vertices, boundary,
                        full_send_dict, ghost_recv_dict, velocity=[1.0, 0.0])

# start the initialisation routines

domain.initialise_visualiser()
#domain.visualise=True
 
#Boundaries

T = Transmissive_boundary(domain)
D = Dirichlet_boundary(array([1.0]))

# ????????????
domain.default_order = 2

# set the boundary information ('exterior; corresponds
# to the ghost triangles sitting along the boundary

domain.set_boundary( {'37': T, 'exterior': T} )

# check the triangularisation

domain.check_integrity()

# set the initial time step

domain.set_quantity('stage', Set_Stage(0.2,0.4,1.0))


# step through time evolving the grid

if myid == 0:
    import time
    t0 = time.time()

domain.visualise = True

for t in domain.evolve(yieldstep = 0.1, finaltime = 2.0):
    if myid == 0:
      domain.write_time()
      
if myid == 0:
    print 'That took %.2f seconds' %(time.time()-t0)
    print 'Communication time %.2f seconds'%domain.communication_time
    print 'Reduction Communication time %.2f seconds'%domain.communication_reduce_time

  
pypar.finalize()