<commit_before>
# ------------------------------------------------------------------------
# Imports
try:
    import fullwv
except ImportError:
    print('Cannot import \'fullwv\'; callback functions will not work!')

import time

# ------------------------------------------------------------------------
# Module variables
starttime = 0

# ------------------------------------------------------------------------
# Callbacks

def CB_the_start():
    global starttime
    fullwv.dbprint('At the start of the program!')
    starttime = time.time()

    fullwv.pythonSolver[:] = True

def CB_the_end():
    fullwv.dbprint('At the end of the program!')
    timediff = time.time() - starttime
    print('\nTime Elapsed: %8.3f\n'% (timediff,))

# ------------------------------------------------------------------------
# Python discretization from Zephyr

import numpy as np
import scipy.sparse
import scipy.sparse.linalg

class Eurus(object):

    c           =   None
    rho         =   None
    nPML        =   None
    freq        =   None
    ky          =   None
    dx          =   None
    dx          =   None
    nx          =   None
    nz          =   None
    freeSurf    =   None

    def __init__(self, systemConfig):

        initMap = {
        #   Argument        Rename to Property
            'c':            None,
            'rho':          None,
            'nPML':         None,
            'freq':         None,
            'ky':           None,
            'dx':           None,
            'dz':           None,
            'nx':           None,
            'nz':           None,
            'freeSurf':     None,
            'theta':        '_theta',
            'eps':          '_eps',
            'delta':        '_delta',
        }

        for key in initMap.keys():
            if key in systemConfig:
                if initMap[key] is None:
                    setattr(self, key, systemConfig[key])
                else:
                    setattr(self, initMap[key], systemConfig[key])

    def _initHelmholtzNinePoint(self):
        """
        An attempt to reproduce the finite-difference stencil and the
        general behaviour of OMEGA by Pratt et al. The stencil is a 9-point
        second-order version based on work by a number of people in the mid-90s
        including Ivan Stekl. The boundary conditions are based on the PML
        implementation by <NAME> in fdfdpml.f.

        SMH 2015: I have modified this code to instead follow the 9-point
        anisotropic stencil suggested by Operto et al. (2009)

        """

        nx = self.nx
        nz = self.nz
        dims = (nz, nx)
        nrows = nx*nz

        c = self.c
        rho = self.rho

        # fast --> slow is x --> y --> z as Fortran

        # Set up physical properties in matrices with padding
        omega   = 2*np.pi * self.freq
        cPad    = np.pad(c, pad_width=1, mode='edge')
        rhoPad  = np.pad(rho, pad_width=1, mode='edge')

        # Horizontal, vertical and diagonal geometry terms
        dx  = self.dx
        dz  = self.dx
        dxx = dx**2
        dzz = dz**2
        dxz = (dxx+dzz)/2
        dd  = np.sqrt(dxz)
        iom = 1j * omega

        # PML decay terms
        # NB: Arrays are padded later, but 'c' in these lines
        #     comes from the original (un-padded) version

        nPML    = self.nPML

        pmldx   = dx*(nPML - 1)
        pmldz   = dz*(nPML - 1)
        pmlr    = 1e-3
        pmlfx   = 3.0 * np.log(1/pmlr)/(2*pmldx**3)
        pmlfz   = 3.0 * np.log(1/pmlr)/(2*pmldz**3)

        dpmlx   = np.zeros(dims, dtype=np.complex128)
        dpmlz   = np.zeros(dims, dtype=np.complex128)
        isnx    = np.zeros(dims, dtype=np.float64)
        isnz    = np.zeros(dims, dtype=np.float64)

        # Only enable PML if the free surface isn't set

        freeSurf = self.freeSurf

        if not freeSurf[2]:
            isnz[-nPML:,:] = -1 # Top

        if not freeSurf[1]:
            isnx[:,-nPML:] = -1 # Right Side

        if not freeSurf[0]:
            isnz[:nPML,:] = 1 # Bottom

        if not freeSurf[3]:
            isnx[:,:nPML] = 1 # Left side

        dpmlx[:,:nPML] = (np.arange(nPML, 0, -1)*dx).reshape((1,nPML))
        dpmlx[:,-nPML:] = (np.arange(1, nPML+1, 1)*dx).reshape((1,nPML))
        dnx     = pmlfx*c*dpmlx**2
        ddnx    = 2*pmlfx*c*dpmlx
        denx    = dnx + iom
        r1x     = iom / denx
        r1xsq   = r1x**2
        r2x     = isnx*r1xsq*ddnx/denx

        dpmlz[:nPML,:] = (np.arange(nPML, 0, -1)*dz).reshape((nPML,1))
        dpmlz[-nPML:,:] = (np.arange(1, nPML+1, 1)*dz).reshape((nPML,1))
        dnz     = pmlfz*c*dpmlz**2
        ddnz    = 2*pmlfz*c*dpmlz
        denz    = dnz + iom
        r1z     = iom / denz
        r1zsq   = r1z**2
        r2z     = isnz*r1zsq*ddnz/denz

        #For now, assume r1x and r1z are the same as xi_x and xi_z from operto et al. (2009)

        #pad edges
        r1x     = np.pad(r1x, pad_width=1, mode='edge')
        r1z     = np.pad(r1z, pad_width=1, mode='edge')

        Xi_x     = 1. / r1x[1,:].reshape((1,nx+2))
        Xi_z     = 1. / r1z[:,1].reshape((nz+2,1))

        # Visual key for finite-difference terms
        # (per Pratt and Worthington, 1990)
        #
        #   This         Original
        # AA BB CC  vs.  AD DD CD
        # DD EE FF  vs.  AA BE CC
        # GG HH II  vs.  AF FF CF

        # Set of keys to index the dictionaries

        # Anisotropic Stencil is 4 times the size, so we define 4 quadrants
        #
        # A =  M1 M2
        #      M3 M4

        M1_keys = ['M1_GG', 'M1_HH', 'M1_II', 'M1_DD', 'M1_EE', 'M1_FF', 'M1_AA', 'M1_BB', 'M1_CC']
        M2_keys = ['M2_GG', 'M2_HH', 'M2_II', 'M2_DD', 'M2_EE', 'M2_FF', 'M2_AA', 'M2_BB', 'M2_CC']
        M3_keys = ['M3_GG', 'M3_HH', 'M3_II', 'M3_DD', 'M3_EE', 'M3_FF', 'M3_AA', 'M3_BB', 'M3_CC']
        M4_keys = ['M4_GG', 'M4_HH', 'M4_II', 'M4_DD', 'M4_EE', 'M4_FF', 'M4_AA', 'M4_BB', 'M4_CC']

        # Diagonal offsets for the sparse matrix formation

        offsets = {
                'GG':   (-1) * nx + (-1),
                'HH':   (-1) * nx + ( 0),
                'II':   (-1) * nx + (+1),
                'DD':   ( 0) * nx + (-1),
                'EE':   ( 0) * nx + ( 0),
                'FF':   ( 0) * nx + (+1),
                'AA':   (+1) * nx + (-1),
                'BB':   (+1) * nx + ( 0),
                'CC':   (+1) * nx + (+1),
        }

        # Need to initialize the PML values

        Xi_x1 = Xi_x[:,0:-2] #left
        Xi_x2 = Xi_x[:,1:-1] #middle
        Xi_x3 = Xi_x[:,2:  ]   #right

        Xi_z1= Xi_z[0:-2,:] #left
        Xi_z2= Xi_z[1:-1,:] #middle
        Xi_z3= Xi_z[2:  ,:] #right<commit_msg>Remove some old lines that don't do anything anymore.<commit_after>
# ------------------------------------------------------------------------
# Imports
try:
    import fullwv
except ImportError:
    print('Cannot import \'fullwv\'; callback functions will not work!')

import time

# ------------------------------------------------------------------------
# Module variables
starttime = 0

# ------------------------------------------------------------------------
# Callbacks

def CB_the_start():
    global starttime
    fullwv.dbprint('At the start of the program!')
    starttime = time.time()

    fullwv.pythonSolver[:] = True

def CB_the_end():
    fullwv.dbprint('At the end of the program!')
    timediff = time.time() - starttime
    print('\nTime Elapsed: %8.3f\n'% (timediff,))

# ------------------------------------------------------------------------
# Python discretization from Zephyr

import numpy as np
import scipy.sparse
import scipy.sparse.linalg

class Eurus(object):

    c           =   None
    rho         =   None
    nPML        =   None
    freq        =   None
    ky          =   None
    dx          =   None
    dx          =   None
    nx          =   None
    nz          =   None
    freeSurf    =   None

    def __init__(self, systemConfig):

        initMap = {
        #   Argument        Rename to Property
            'c':            None,
            'rho':          None,
            'nPML':         None,
            'freq':         None,
            'ky':           None,
            'dx':           None,
            'dz':           None,
            'nx':           None,
            'nz':           None,
            'freeSurf':     None,
            'theta':        '_theta',
            'eps':          '_eps',
            'delta':        '_delta',
        }

        for key in initMap.keys():
            if key in systemConfig:
                if initMap[key] is None:
                    setattr(self, key, systemConfig[key])
                else:
                    setattr(self, initMap[key], systemConfig[key])

    def _initHelmholtzNinePoint(self):
        """
        An attempt to reproduce the finite-difference stencil and the
        general behaviour of OMEGA by Pratt et al. The stencil is a 9-point
        second-order version based on work by a number of people in the mid-90s
        including <NAME>. The boundary conditions are based on the PML
        implementation by <NAME> in fdfdpml.f.

        SMH 2015: I have modified this code to instead follow the 9-point
        anisotropic stencil suggested by Operto et al. (2009)

        """

        nx = self.nx
        nz = self.nz
        dims = (nz, nx)
        nrows = nx*nz

        c = self.c
        rho = self.rho

        # fast --> slow is x --> y --> z as Fortran

        # Set up physical properties in matrices with padding
        omega   = 2*np.pi * self.freq
        cPad    = np.pad(c, pad_width=1, mode='edge')
        rhoPad  = np.pad(rho, pad_width=1, mode='edge')

        # Horizontal, vertical and diagonal geometry terms
        dx  = self.dx
        dz  = self.dx
        dxx = dx**2
        dzz = dz**2
        dxz = (dxx+dzz)/2
        dd  = np.sqrt(dxz)
        iom = 1j * omega

        # PML decay terms
        # NB: Arrays are padded later, but 'c' in these lines
        #     comes from the original (un-padded) version

        nPML    = self.nPML

        pmldx   = dx*(nPML - 1)
        pmldz   = dz*(nPML - 1)
        pmlr    = 1e-3
        pmlfx   = 3.0 * np.log(1/pmlr)/(2*pmldx**3)
        pmlfz   = 3.0 * np.log(1/pmlr)/(2*pmldz**3)

        dpmlx   = np.zeros(dims, dtype=np.complex128)
        dpmlz   = np.zeros(dims, dtype=np.complex128)
        isnx    = np.zeros(dims, dtype=np.float64)
        isnz    = np.zeros(dims, dtype=np.float64)

        # Only enable PML if the free surface isn't set

        freeSurf = self.freeSurf

        if not freeSurf[2]:
            isnz[-nPML:,:] = -1 # Top

        if not freeSurf[1]:
            isnx[:,-nPML:] = -1 # Right Side

        if not freeSurf[0]:
            isnz[:nPML,:] = 1 # Bottom

        if not freeSurf[3]:
            isnx[:,:nPML] = 1 # Left side

        dpmlx[:,:nPML] = (np.arange(nPML, 0, -1)*dx).reshape((1,nPML))
        dpmlx[:,-nPML:] = (np.arange(1, nPML+1, 1)*dx).reshape((1,nPML))
        dnx     = pmlfx*c*dpmlx**2
        ddnx    = 2*pmlfx*c*dpmlx
        denx    = dnx + iom
        r1x     = iom / denx
        r1xsq   = r1x**2
        r2x     = isnx*r1xsq*ddnx/denx

        dpmlz[:nPML,:] = (np.arange(nPML, 0, -1)*dz).reshape((nPML,1))
        dpmlz[-nPML:,:] = (np.arange(1, nPML+1, 1)*dz).reshape((nPML,1))
        dnz     = pmlfz*c*dpmlz**2
        ddnz    = 2*pmlfz*c*dpmlz
        denz    = dnz + iom
        r1z     = iom / denz
        r1zsq   = r1z**2
        r2z     = isnz*r1zsq*ddnz/denz

        #For now, assume r1x and r1z are the same as xi_x and xi_z from operto et al. (2009)

        #pad edges
        r1x     = np.pad(r1x, pad_width=1, mode='edge')
        r1z     = np.pad(r1z, pad_width=1, mode='edge')

        Xi_x     = 1. / r1x[1,:].reshape((1,nx+2))
        Xi_z     = 1. / r1z[:,1].reshape((nz+2,1))

        # Visual key for finite-difference terms
        # (per Pratt and Worthington, 1990)
        #
        #   This         Original
        # AA BB CC  vs.  AD DD CD
        # DD EE FF  vs.  AA BE CC
        # GG HH II  vs.  AF FF CF

        # Set of keys to index the dictionaries

        # Anisotropic Stencil is 4 times the size, so we define 4 quadrants
        #
        # A =  M1 M2
        #      M3 M4

        # Diagonal offsets for the sparse matrix formation

        offsets = {
                'GG':   (-1) * nx + (-1),
                'HH':   (-1) * nx + ( 0),
                'II':   (-1) * nx + (+1),
                'DD':   ( 0) * nx + (-1),
                'EE':   ( 0) * nx + ( 0),
                'FF':   ( 0) * nx + (+1),
                'AA':   (+1) * nx + (-1),
                'BB':   (+1) * nx + ( 0),
                'CC':   (+1) * nx + (+1),
        }

        # Need to initialize the PML values

        Xi_x1 = Xi_x[:,0:-2] #left
        Xi_x2 = Xi_x[:,1:-1] #middle
        Xi_x3 = Xi_x[:,2:  ]   #right

        Xi_z1= Xi_z[0:-2,:] #left
        Xi_z2= Xi_z[1:-1,:] #middle
        Xi_z3= Xi_z[2:  ,:] #right
{'fraction_non_alphanumeric': 0.16864881454768665, 'fraction_numerical': 0.02608762372439193, 'mean_word_length': 1.3003882809742322, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3886729', 'n_tokens_mistral': 4809, 'n_tokens_neox': 4494, 'n_words': 1164}