mmdgAnalysis.py

import numpy as np
from os import mkdir, remove, system
from os.path import join, exists
from subprocess import call
from math import log


# changes the parameter file such that gridfile is used
def writeINI(gridfile):
    # storage for file content
    lines = None

    # new line to write
    newline = "gridfile = " + gridfile

    # indicates whether the parameter gridfile already exists
    hasGridlineParameter = False

    # read file content and replace gridfile parameter if existent
    with open("parameterMMDG.ini", "r") as file:
        lines = file.read().splitlines()

        for i in range(len(lines)):
            if (lines[i].startswith("gridfile")):
                lines[i] = newline
                hasGridlineParameter = True
                break

    # add parameter gridfile if not already existent
    if (not hasGridlineParameter):
        lines.append(newline)

    # write parameters to file
    with open("parameterMMDG.ini", "w") as file:
        file.write("\n".join(lines))


# reads error, run time and maximum edge length, i.e. the output from the
# excuted discontinuous Galerkin scheme
def readData():
    file = open("convergenceData", "r")

    # read data
    errorBulk = float(file.readline())
    errorInterface = float(file.readline())
    errorTotal = float(file.readline())
    hMax = float(file.readline())
    hMaxInterface = float(file.readline())
    timeTotal = float(file.readline())

    file.close()

    # return data
    return np.array([errorBulk, errorInterface, errorTotal, hMax, \
        hMaxInterface, timeTotal])


# determines the problem type and the parameter "fastEOC" from the ini file
def getProblemType ():
    # storage for problem type (1 to 6)
    problemType = None

    # boolean that determines whether to exclude computationally expensive grids
    fastEOC = False

    # read file content and get problem type
    with open("parameterMMDG.ini", "r") as file:
        lines = file.read().splitlines()

        for i in range(len(lines)):
            if (lines[i].startswith("problem")):
                problemType = int(lines[i].rstrip()[-1])
            elif (lines[i].startswith("fastEOC") and \
             lines[i].rstrip().lower().endswith("true")):
                fastEOC = True

    return problemType, fastEOC


# determines experimental order of convergence (EOC),
# data is assumed to be a numpy array with
# first column: error,
# second column: maximum edge length
def EOCloop (data):
    EOC = np.empty(data.shape[0] - 1)

    for i in range(len(EOC)):
        if data[i, 0] != 0:
            EOC[i] = log(data[i, 0] / data[i+1, 0]) / \
                log(data[i, 1] / data[i+1, 1])
        else:
            EOC[i] = -1.0 * float('inf')
    # end for i

    return EOC


# === main ===

# create directory for the plots if necessary
if not exists("plots"):
    mkdir("plots")

# number of runs from which the run times are taken as mean value
numberOfMeans = 1

# names of the grid files sorted by refinement
# 2d characteristic lengths: np.reciprocal(np.logspace(0,2.7,12))
# 3d characteristic lengths: np.reciprocal(np.logspace(0,1.55,7))
gridfiles_2d = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L"]
gridfiles_3d = ["A", "B", "C", "D", "E", "F", "G"]

# get problem type (integer from 1 to 6) from ini file
problemType, fastEOC = getProblemType()

if problemType in [1, 4, 6, 7]: # these problems use the same grids
    problemType = 2

if fastEOC: # exclude computationally expensive grids
    gridfiles_2d = gridfiles_2d[:-3]
    gridfiles_3d = gridfiles_3d[:-2]
else:
    print("fastEOC = false. This may take some minutes.\n")

gridfiles_2d = ["mmdg" + str(problemType) + "_" + s for s in gridfiles_2d]
gridfiles_3d = ["mmdg" + str(problemType) + "_" + s for s in gridfiles_3d]

for dim in [2, 3]:

    gridfiles = gridfiles_2d if (dim == 2) else gridfiles_3d

    # storage for errors, run times and largestEdgeLength
    # column 0: bulk error
    # column 1: interface error
    # column 2: total error
    # column 3: maximum edge length (bulk)
    # column 4: maximum edge length (interface)
    # column 5: total run time
    # different rows correspond to different numbers of grid elements
    data = np.zeros((len(gridfiles), 6))

    for i in range(len(gridfiles)):
        # write dgf file
        writeINI(gridfiles[i])

        # take mean value over numberOfMeans runs
        for j in range(numberOfMeans):
            # run executable
            call(join(".", "mmdg-" + str(dim) + "d"))

            # store data
            data[i, :] += readData()
        # end for j
    # end for i

    # divide by numberOfMeans
    data /= numberOfMeans

    # determine experimental order of convergence (EOC)
    EOCbulk = EOCloop(data[:, [0,3]])
    EOCinterface = EOCloop(data[:, [1,4]])
    EOCtotal = \
        EOCloop(np.column_stack([data[:,2], np.maximum(data[:,3], data[:,4])]))

    # print EOCs
    print("\n-------------------------------------------------------------\n")
    print("EOCs (" + str(dim) + "d):")
    print("\nbulk:")
    print(EOCbulk)
    print("\ninterface:")
    print(EOCinterface)
    print("\ntotal:")
    print(EOCtotal)
    print("\n=============================================================\n\n")

    # write data to file
    np.savetxt(join("plots", "mmdgAnalysis_" + str(dim) + "d"), data)
    np.savetxt(join("plots", "mmdgEOC_" + str(dim) + "d"), EOCtotal)
# end for dim

# remove output files
remove("convergenceData")