convergencetest.py 2.86 KB
Newer Older
1
#!/usr/bin/env python3
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

from math import *
import subprocess
import sys

if len(sys.argv) < 2:
    sys.stderr.write('Please provide a single argument <testname> to the script\n')
    sys.exit(1)

testname = str(sys.argv[1])
testargs = [str(i) for i in sys.argv][2:]

# remove the old log file
subprocess.call(['rm', testname + '.log'])
print("Removed old log file ({})!".format(testname + '.log'))

# do the runs with different refinement
for i in [0, 1]:
    subprocess.call(['./' + testname] + testargs + ['-Grid.Refinement', str(i),
                                      '-Problem.Name', testname])

# check the rates and append them to the log file
logfile = open(testname + '.log', "r+")

errorP = []
errorVx = []
errorVy = []
for line in logfile:
    line = line.strip("\n")
    line = line.strip("\[ConvergenceTest\]")
    line = line.split()
    errorP.append(float(line[2]))
    errorVx.append(float(line[5]))
    errorVy.append(float(line[8]))

resultsP = []
resultsVx = []
resultsVy = []
logfile.truncate(0)
logfile.write("n\terrorP\t\trateP\t\terrorVx\t\trateVx\t\terrorVy\t\trateVy\n")
logfile.write("-"*50 + "\n")
for i in range(len(errorP)-1):
    if isnan(errorP[i]) or isinf(errorP[i]):
        continue
    if not ((errorP[i] < 1e-12 or errorP[i+1] < 1e-12) and (errorVx[i] < 1e-12 or errorVx[i+1] < 1e-12) and (errorVy[i] < 1e-12 or errorVy[i+1] < 1e-12)):
        rateP = (log(errorP[i])-log(errorP[i+1]))/log(2)
        rateVx = (log(errorVx[i])-log(errorVx[i+1]))/log(2)
        rateVy = (log(errorVy[i])-log(errorVy[i+1]))/log(2)
        message = "{}\t{:0.4e}\t{:0.4e}\t{:0.4e}\t{:0.4e}\t{:0.4e}\t{:0.4e}\n".format(i, errorP[i], rateP,  errorVx[i], rateVx, errorVy[i], rateVy)
        logfile.write(message)
        resultsP.append(rateP)
        resultsVx.append(rateVx)
        resultsVy.append(rateVy)
    else:
        logfile.write("error: exact solution!?")
i = len(errorP)-1
message = "{}\t{:0.4e}\t\t{}\t{:0.4e}\t\t{}\t{:0.4e}\t\t{}\n".format(i, errorP[i], "",  errorVx[i], "", errorVy[i], "")
logfile.write(message)

logfile.close()
print("\nComputed the following convergence rates for {}:\n".format(testname))
subprocess.call(['cat', testname + '.log'])

Timo Koch's avatar
Timo Koch committed
65
66
67
def mean(numbers):
    return float(sum(numbers)) / len(numbers)

68
# check the rates, we expect rates around 2
Timo Koch's avatar
Timo Koch committed
69
if mean(resultsP) < 2.05 and mean(resultsP) < 1.95:
70
71
72
    sys.stderr.write("*"*70 + "\n" + "The convergence rates for pressure were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
    sys.exit(1)

Timo Koch's avatar
Timo Koch committed
73
if mean(resultsVx) < 2.05 and mean(resultsVx) < 1.95:
74
75
76
    sys.stderr.write("*"*70 + "\n" + "The convergence rates for x-velocity were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
    sys.exit(1)

Timo Koch's avatar
Timo Koch committed
77
if mean(resultsVy) < 2.05 and mean(resultsVy) < 1.95:
78
79
80
81
    sys.stderr.write("*"*70 + "\n" + "The convergence rates for y-velocity were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
    sys.exit(1)

sys.exit(0)