vergleich.py

Sfdtd Simple finite-difference time-domain

# vergleich.py
#
# Some comparison of vergleich*.plt, verlauf*.plt, cube*.data 
# and poynting*.data output files from Sfdtd
#
# (This is a quite dirty python code... )
#
#    Copyright (C) 2007  Paul Panserrieu, < [e-mail] >
#
#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License.
#
# last modified: 26-10-2007 05:59:54 PM CEST

import os

try:
    import Numeric
except ImportError:
    import sys
    # if you have installed Numeric in a different directory, change the path
    sys.path.append('/home/p/peutetre/bin/lib/python/Numeric')
    try:
        import Numeric
    except ImportError:
        print "Numeric module not available..."
        pass

print """ >~~~~~~~~~~~~~~~~>>~~~~~~~~>>>~~~~~~~~>>>> """
print """ python vergleich.py """
print """ Vergleich von Dateien aus Sfdtd """
print """ Comparison of files from Sfdtd """
print """ Paul Panserrieu <[e-mail]> """
print """ >~~~~~~~~~~~~~~~~>>~~~~~~~~>>>~~~~~~~~>>>> """

file_id_1 = '_1_vergleich_%s_00%i0%i.plt'
file_id_2 = '_2_vergleich_%s_00%i0%i.plt'
verlauf_id_1 = '_1_verlauf2d_%s_0090%i.plt'
rela_id_1 = '_1_rela_%s_0090%i.plt'

def w():
    print os.getcwd()

def l():
    print os.listdir('.')

def cd(path='~'):
    os.chdir(path)

def up():
    os.chdir('..')

def down():
    files = os.listdir('.') 
    dirs = []
    for i in files:
      if os.path.isdir(i):
          dirs.append(i)
    n = 0
    for i in dirs:
        print n, i
        n += 1
    if n > 0:
        print 'change dir to ?'
        n = raw_input(':  ')
        if n in range(0, len(dirs)):
            os.chdir(os.getcwd() + '/' + str(dirs[n]))

class curve:

    def __init__(self):
        self.data = dict()
        self.scale = dict()
        self.title = ""
        return None

    def __repr__(self):
        chaine = "keys: \n"
        for i in self.data.keys():
            if i != self.data.keys()[-1]: 
                chaine += ' %s\n' %(i)
            else:
                chaine += ' %s' %(i)
        return chaine

    def show(self, key):
        return self.data[key]

    def set_scale(self, key, factor):
        if self.data.has_key(key):
            self.scale[key] = factor
            return None
        else:
            return None
  
    def set_title(self, title):
        self.title = title

    def write(self, key_list, file_name, line_title=""):
        for i in key_list:
            if not self.data.has_key(i):
                print 'key %s does not exist' %i
                return None
        try:
            f = open(file_name, 'w')
        except IOError:
            print "does not manage to create file %s" % file_name
            return None
        
        f.write("$ DATA = CURVE2D \n")
        f.write("%"+" toplabel= '%s' \n" % self.title)
        f.write("% xlabel= 'time-step'\n")
        f.write("% ylabel= ' ' \n")
        a = 1

        for i in key_list:
            if line_title != "" and int(i.split('_')[1]) == 1:
                f.write("%"+" linelabel = %s \n" % line_title)
            else:
                try:
                    f.write("%"+" linelabel = %s \n" % i.split('_')[3])
                except IndexError:
                    f.write("%"+" linelabel = %s \n" % i)
            f.write("%"+" linecolor = %i linetype = 1 \n" % a)
            b = 0; scale = 1.0
            if i in self.scale.keys():
                scale = self.scale[i]
            for j in self.data[i]:
                f.write("%i  %s\n" %(b*scale, j))
                b += 1
            f.write("\n")
            a += 1
        f.write(" $END \n")
        f.close()

    def read(self, file_name):
        try:
            f = open(file_name, 'r')
        except IOError:
            print "%s cannot be opened..." % file_name
            return 13
        chaine = ""
        self.data['_1_' + file_name] = []
        for i in range(1,9):
             chaine = f.next()
        while not chaine.startswith('\n') or chaine.find('$') == -1:
            try:
                try:
                    self.data['_1_' + file_name].append(float(chaine.split()[-1]))
                    chaine = f.next()
                except ValueError:
                    break
            except IndexError:
                break
        if chaine.startswith(' $ END'):
            return None
        else:
            self.data['_2_' + file_name] = []
            for i in range(1,4):
                chaine = f.next()
            while chaine.find('$') == -1:
                try:
                    self.data['_2_' + file_name].append(float(chaine.split()[-1]))
                    chaine = f.next()
                except IndexError:
                    break
        f.close()
        return None

    def add(self, array_id, array):
        self.data[array_id]=array

    def find_max(self, array_id, offset = 0):
        if array_id in self.data.keys():
            pass
        else: return None
        tmp = 0.0e0
        for i in self.data[array_id][offset:]:
            if abs(float(i)) > abs(tmp):
                tmp = float(i)
        return tmp  

    def subtract(self, first_id, second_id, new_id):

        if (len(self.data[first_id]) == len(self.data[second_id])):
            pass
        else:
            print '%s and %s do not have the same size' %(first_id, second_id)
            return None

        tmp = []
        first = self.show(first_id)
        second = self.show(second_id)
        
        for i in range(0, len(first)):
            tmp.append(first[i]-second[i])

        self.add(new_id, tmp)
        return None

    def sub_rela(self, first_id, second_id, new_id, maxi):

        if (len(self.data[first_id]) == len(self.data[second_id])):
            pass
        else:
            print '%s and %s do not have the same size' %(first_id, second_id)
            return None

        tmp = []
        first = self.show(first_id)
        second = self.show(second_id)

        for i in range(0, len(first)):
            tmp.append(100.0*abs(first[i]-second[i])/maxi)
        self.add(new_id, tmp)
        return None 

    def norm(self, d, new_id):
        tmp = []
        first = self.show(d)
        maxi = self.find_max(d)
        for i in first:
            tmp.append(i/maxi)
        self.add(new_id, tmp)
        

def error_1(offset=0):

    """ local error: numerical result v.s. analytical solution """

    a = os.listdir('.')
    files = []
    for i in a:
        if i.endswith('.plt') and i.startswith('ver'):
            files.append(i)
    del(a)
    ids = []
    for i in files:
        tmp = i.split('_')[1]
        if tmp not in ids:
            ids.append(tmp)
    f = curve()
    maxi = dict()
    for i in files:
        f.read(i)
    for j in [1,2]:
        for k in [1,2,3]:
            plots = []
            for i in ids:
                if ids.index(i) == 0:
                    plots.append(file_id_1%(i,j,k))
                    maxi[file_id_1%(i,j,k)]=\
                    (f.find_max(file_id_1%(i,j,k), offset))
                    plots.append(file_id_2%(i,j,k))
                else:
                    maxi[file_id_1%(i,j,k)]=\
                    (f.find_max(file_id_1%(i,j,k), offset))
                    plots.append(file_id_2%(i,j,k))
            f.write(plots, 'all_%i0%i.plt'%(j,k), 'analytic')         
    for k in range(1,6):
        plots = []
        for i in ids:
            tmp = f.show(verlauf_id_1%(i,k))
            tmp1 = []
            for j in tmp:
                tmp1.append(100.0*abs(j)\
                            /abs(maxi[file_id_1%(i,oups(k),chut(k))]))
            f.add(rela_id_1%(i,k), tmp1)
            plots.append(rela_id_1%(i,k))
            if k == 3:
                print 'abso', i, f.find_max(verlauf_id_1%(i,k),offset)
                print 'rela', i, f.find_max(rela_id_1%(i,k), offset)
        f.write(plots, 'all_error_rela_analytic_1_90%i.plt'%(k))
    return None

def error_1_bis(ref_id, abso = 1):

    """ local error: numerical result v.s. reference solution """

    a = os.listdir('.')
    files = []
    for i in a:
        if i.endswith('.plt') and i.startswith('vergleich'):
            files.append(i)
    del(a)
    ids = []
    for i in files:
        tmp = i.split('_')[1]
        if tmp not in ids:
            ids.append(tmp)
    if ref_id not in ids:
        print 'ref_id:%s is not available' % ref_id
        return 13
    f = curve()
    for i in files:
        f.read(i)
    for j in [1,2]:
        for k in [1,2,3]:
            plots = []
            for i in ids:
                if i == ref_id:
                    if abso == 1:
                        pass
                    else:
                        maxi = 0.0
                        for l in f.show(file_id_2%(ref_id,j,k)):
                            if abs(l) > maxi: maxi = abs(l)
                else: pass
            for i in ids:
                if i !=  ref_id:
                    if abso == 1:
                        f
                        f.subtract(  file_id_2%(i,j,k)\
                                   , file_id_2%(ref_id,j,k)\
                                   , 'sub_%s_00%i0%i.plt'%(i,j,k) )  
                        plots.append('sub_%s_00%i0%i.plt'%(i,j,k))
                    else:
                        f.sub_rela(  file_id_2%(i,j,k)\
                                   , file_id_2%(ref_id,j,k)\
                                   , 'sub_%s_00%i0%i.plt'%(i,j,k)\
                                   , maxi )
                        plots.append('sub_%s_00%i0%i.plt'%(i,j,k))
                else:
                    pass
            if abso == 1:
                f.write(plots, 'sub_abs_%i0%i.plt'%(j,k))
            else:
                f.write(plots, 'sub_rela_%i0%i.plt'%(j,k))
    return None
            
class cube:

    def __init__(self):
        self.limit = 0
        return None

    def read(self, file_name):
        try:
            f = open(file_name, 'r')
        except IOError:
            print "%s cannot be opened..." % file_name
            return 13
        self.limit = int(f.next())
        self.data = Numeric.zeros((2*self.limit+1 \
                                 , 2*self.limit+1 \
                                 , 2*self.limit)  \
                                 , typecode='fd')
        f.next()
        for k in range(-self.limit, self.limit):
            for j in range(-self.limit, self.limit+1):
                for i in range(-self.limit, self.limit+1):
                    self.data[i][j][k] = float(f.next().split()[-1])
        f.close()
        return None

def error_2(ref_id):

    """ global error: numerical result v.s. reference solution """

    a = os.listdir('.')
    files = []
    for i in a:
        if i.endswith('.data') and i.startswith('cube'):
            files.append(i)
    del(a)
    files.sort()
    ids = []
    for i in files:
        tmp = i.split('_')[1]
        if not tmp in ids:
            ids.append(tmp)
    if ref_id not in ids:
        print 'ref_id:%s is not available' % ref_id
        return 13
    maxi = dict()
    for i in ids:
        maximum = 0
        for j in files:
            if j.split('_')[1] == i:
                if int(j.split('_')[-1].split('.')[0]) > maximum:
                    maximum = int(j.split('_')[-1].split('.')[0])
        maxi[i] = maximum
    max_timestep = min(maxi.values())
    ids.remove(ref_id)
    data = dict()
    print "available ids: %s" %ids
    max_ref=[]
    for t in range(1, max_timestep+1):
        ref = cube()
        ref.read('cube_%s_%s.data' %(ref_id, str(t).zfill(5)))
        tmp1 = 0.0
        for i in range(-ref.limit, ref.limit+1):
                for j in range(-ref.limit, ref.limit+1):
                    for k in range(-ref.limit, ref.limit):
                        tmp1 += pow(ref.data[i][j][k],2)
        max_ref.append(tmp1)
        for d in ids:
            tmp = cube()
            tmp.read('cube_%s_%s.data' %(d, str(t).zfill(5)))
            value = 0.0
            for i in range(-ref.limit, ref.limit+1):
                for j in range(-ref.limit, ref.limit+1):
                    for k in range(-ref.limit, ref.limit):
                        value += pow(tmp.data[i][j][k] - ref.data[i][j][k],2)
            del(tmp)
            if t == 1:
                data[d] = [0.0]
            data[d].append(value)
        del(ref)
        print "time step %i done" %t
    m_ref = max(max_ref)
    for d in ids:
        tmp = []
        for n in data[d]:
            tmp.append(100.0*n/m_ref)
        data[d] = tmp
    plot = curve()
    for d in data.keys():
        plot.add(d,data[d])
    plot.write(data.keys(), 'global_error.plt')
    del(plot)
    return None

def chut(a):
  if a < 4: return a
  else: return a-3

def oups(a):
  if a in [1,2,3]: return 1
  else: return 2

class fluss:

    def __init__(self):
        self.data = []
        return None

    def read(self, file_name):
        try:
            f = open(file_name, 'r')
        except IOError:
            print "%s cannot be opened..." % file_name
            return 13
        f.next()
        tmp = f.next()
        self.data.append(float(tmp.split()[-1]))
        f.close()
        return None


def error_3(ref_id):

    """ energy flow error: numerical result v.s. reference solution """

    a = os.listdir('.')
    files = []
    for i in a:
        if i.endswith('.data') and i.startswith('poynting'):
            files.append(i)
    del(a)
    files.sort()
    ids = []
    for i in files:
        tmp = i.split('_')[1]
        if not tmp in ids:
            ids.append(tmp)
    if ref_id not in ids:
        print 'ref_id:%s is not available' % ref_id
        return 13
    maxi = dict()
    for i in ids:
        maximum = 0
        for j in files:
            if j.split('_')[1] == i:
                if int(j.split('_')[-1].split('.')[0]) > maximum:
                    maximum = int(j.split('_')[-1].split('.')[0])
        maxi[i] = maximum
    max_timestep = min(maxi.values())
    ids.remove(ref_id)
    data = dict()
    print "available ids: %s" %ids
    max_ref= 0.0
    ref = fluss()
    ref_data = []
    result = dict()
    for t in range(0, max_timestep):
        ref.read('poynting_%s_%s.data' %(ref_id, str(t).zfill(5)))
        if abs(ref.data[t]) > abs(max_ref):
            max_ref = ref.data[t]
    print ref_id, 'max abs(flow):', max_ref

    ref_data = ref.data
    del(ref)
    for d in ids:
        tmp = fluss()
        for t in range(0, max_timestep):
            tmp.read('poynting_%s_%s.data' %(d, str(t).zfill(5)))
        data[d] = tmp.data 
        del(tmp)

    for d in ids:
        tmp = []
        for t in range(0, max_timestep):
            tmp.append(100.0*abs(data[d][t] - ref_data[t])/abs(max_ref))
        result[d] = tmp
        print d, max(result[d])
    plot = curve()
    for d in result.keys():
        plot.add(d,result[d])
    plot.write(result.keys(), 'energy_error.plt')
    del(plot)
    return None