source: trunk/anuga_work/anuga_cuda/src/anuga_HMPP/hmpp_domain.py @ 9329

# -*- coding: utf-8 -*-


import ctypes
import sys


import numpy


# This is used to fix the dlopen() issues for python
# and make the reference links (like OpenHMPP global)
flags = sys.getdlopenflags()
sys.setdlopenflags(flags | ctypes.RTLD_GLOBAL)

from hmpp_python_glue import *
sys.setdlopenflags(flags)


from anuga import Domain
from anuga import Reflective_boundary
from anuga import Dirichlet_boundary
from anuga import Time_boundary
from anuga import Transmissive_boundary

import numpy as num



class HMPP_domain(Domain):
    """ANUGA in OpenHMPP implementation.
    
    In this version, only the mesh information is generated by Python 
    grogram, while the overall evolve procedure is moved to C (OpenHMPP)
    implementation. This is for the benefits of the performance, since 
    passing host control back and force between Python and C throught 
    Python/C API can be costly.
    """

    def __init__(self, 
                coordinates=None, 
                vertices=None,
                 boundary=None,
                 source=None,
                 triangular=None,
                 conserved_quantities=None,
                 evolved_quantities=None,
                 other_quantities=None,
                 tagged_elements=None,
                 geo_reference=None,
                 use_inscribed_circle=False,
                 mesh_fulename=None,
                 use_cache=False,
                 verbose=False,
                 full_send_dict=None,
                 ghost_recv_dict=None,
                 starttime=0.0,
                 processor=0,
                 numproc=1,
                 number_of_full_nodes=None,
                 number_of_full_triangles=None,
                 ghost_layer_width=2
                 ):
        """The init routain of the ANUGA to create the mesh."""

        Domain.__init__(self,
                        coordinates,
                        vertices,
                        boundary,
                        full_send_dict=full_send_dict,
                        ghost_recv_dict=ghost_recv_dict,
                        number_of_full_nodes=number_of_full_nodes,
                        number_of_full_triangles=number_of_full_triangles,
                        geo_reference=geo_reference) #jj added this


    def convert_boundary_elements(self):
        """Write down all mesh boundary information into file.
        
        This is used to pass boundary information to C implementation.
        """

        fileHandle = open('boundary_names', 'w')

        boundary_cnt = 0
        for name in self.tag_boundary_cells:
            B = self.boundary_map[name]
            if B is None:
                continue


            fileHandle.write("%s\n" % name)
            if isinstance(B, Reflective_boundary):
                fileHandle.write("0\n")
            elif isinstance(B, Dirichlet_boundary):
                fileHandle.write("1\n")
            else:
                print B

            boundary_cnt += 1

            if name == 'open':
                self.openArr = numpy.asarray(
                        self.tag_boundary_cells[name], dtype=numpy.int64)
            elif name == 'exterior':
                self.exterior = numpy.asarray(
                        self.tag_boundary_cells[name], dtype=numpy.int64)

        fileHandle.close()
        
        return boundary_cnt



    def evolve(self,
                yieldstep=0.0,
                finaltime=1.0,
                duration=0.0,
                skip_initial_step=False):
        """Evolve function
        
        Here the overall evolution starts, but the procedures are done in
        the C implementation. 
        """

        if self.store is True and self.get_time() == self.get_starttime():
            self.initialise_storage()

        from hmpp_python_glue import hmpp_evolve
        from anuga.config import epsilon

        if self.timestepping_method == 'euler':
            timestepping_method = 1
        elif self.timestepping_method == 'rk2':
            timestepping_method = 2
        elif self.timestepping_method == 'rk3':
            timestepping_method = 3
        else:
            timestepping_method = 4


        if self.flow_algorithm == 'tsunami':
            flow_algorithm = 1
        elif self.flow_algorithm == 'yusuke':
            flow_algorithm = 2
        else:
            flow_algorithm = 3
       

        if self.compute_fluxes_method == 'original':
            compute_fluxes_method = 0
        elif self.compute_fluxes_method == 'wb_1':
            compute_fluxes_method = 1
        elif self.compute_fluxes_method == 'wb_2':
            compute_fluxes_method = 2
        elif self.compute_fluxes_method == 'wb_3':
            compute_fluxes_method = 3
        elif self.compute_fluxes_method == 'tsunami':
            compute_fluxes_method = 4
        else:
            compute_fluxes_method = 5

    
        boundary_cnt = self.convert_boundary_elements()                
            
        import time 
        ini_time = time.time()

        yield_step = 0
        while True :
            tmp_timestep = hmpp_evolve(self,
                    yieldstep,
                    finaltime,
                    duration,
                    epsilon,
                    skip_initial_step,

                    numpy.int32( compute_fluxes_method ),
                    numpy.int32( flow_algorithm ),
                    numpy.int32( timestepping_method ),
                    #FIXME: get issues on parsing the boundary 
                    numpy.int64( 0 ),
                    numpy.int32( yield_step )
                    )

            yield_step = 1
            print " Python: tmp_timestep %lf " % tmp_timestep

            if tmp_timestep >= finaltime - epsilon: 
                fin_time = time.time()
                break