# 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, < peutetre@cs.tu-berlin.de >##    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 CESTimport ostry:    import Numericexcept 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..."        passprint """ >~~~~~~~~~~~~~~~~>>~~~~~~~~>>>~~~~~~~~>>>> """print """ python vergleich.py """print """ Vergleich von Dateien aus Sfdtd """print """ Comparison of files from Sfdtd """print """ Paul Panserrieu <peutetre@cs.tu-berlin.de> """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 Nonedef 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 Nonedef 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 Nonedef chut(a):  if a < 4: return a  else: return a-3def oups(a):  if a in [1,2,3]: return 1  else: return 2class 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 Nonedef 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