Skip to content
Snippets Groups Projects

WIP: Add test computation for bending isometries using reducedcubichermitetrianglebasis

Open WIP: Add test computation for bending isometries using reducedcubichermitetrianglebasis
feature/bendingIsometries into master
Open Klaus Böhnlein requested to merge feature/bendingIsometries into master
Viewing commit eb695741
Show latest version
1 file
+ 430
− 0
Compare changes
  • Side-by-side
  • Inline
bending-isometries-testsuite/bending-isometries-testsuite.py 0 → 100644
+ 430
− 0
import os
import subprocess
import sys
import re
import numpy as np
from prettytable import PrettyTable
# import latextable
# from texttable import Texttable
from tabulate import tabulate
# import time
import importlib
####################################################################
# Usage:
#
# python3 harmonicmaps-numerical-tests-batch.py TAURUS_PATH_Flag <slurm_array_task_id>
#
####################################################################
# def computeHarmonicMap(minLevel,maxLevel,dim,targetDim,order,interpolationMethod,maxTrustRegionSteps,randomInitialIterate,readConfiguration,configurationFile,perturbation,perturbationRadius,epsilon,tolerance,executable,parameterFile,pythonPath,resultPath):
def computeBendingIsometry(executable,parameterFile,pythonPath,resultPath,instrumentedPath,gridLevels, conforming_DiscreteJacobian ,INSTRUMENTED):
# print('-----Run computeHarmonicMap----')
processList = []
for level in range(gridLevels[0],gridLevels[1]+1):
# print("Run harmonicmaps_intoR"+ str(parameterSet.targetDim) +" with order "+str(parameterSet.order) +" on grid level " + str(level) )
# print("Run harmonicmaps_intoR"+ str(parameterSet.targetDim) +" with order "+str(order) +" on grid level " + str(level) )
# LOGFILE = "./harmonicmapsR"+ str(dim) +"_intoR"+ str(targetDim) + "_deg" + str(kappa) + "_"+ interpolationMethod + "_" + str(order) + "_" + str(level) + ".log"
# if int(sys.argv[1]):
# LOGFILE = parameterFile + "_order" + str(order) + "_" + interpolationMethod + "_level" + str(level) + ".log"
# else:
if conforming_DiscreteJacobian:
conformity = "_conforming"
else:
conformity = "_nonconforming"
LOGFILE = resultPath + "/" + parameterFile + conformity + "_level" + str(level) + ".log"
print('LOGFILE:',LOGFILE)
print('executable:',executable)
print('parameterFile:', parameterFile)
print('resultPath:', resultPath)
# print('NOW RUN PROGRAM on level:' + str(level))
# start_time = time.time()
p = subprocess.Popen(executable + " " + pythonPath + " " + parameterFile
+ " -numLevels " + str(level)
+ " -resultPath " + str(resultPath)
+ " -instrumentedPath " + str(instrumentedPath)
+ " -conforming_DiscreteJacobian " + str(conforming_DiscreteJacobian)
+ " -instrumented " + str(INSTRUMENTED)
+ " | tee " + LOGFILE, shell=True)
p.wait() # wait
# print("--- %s seconds ---" % (time.time() - start_time))
print('------FINISHED PROGRAM on level:' + str(level))
processList.append(p)
# Wait for all simulation subprocesses before proceeding to the error measurement step
exit_codes = [p.wait() for p in processList]
return
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# # for pretty table output install :
# - python3-ptable
# - python3-prettytable
# 1. experiment, 2. conforming_DiscreteJacobian_, 3. gridLevels, 4. INSTRUMENTED
scenarios = [["L-clamped-Plate" , 1, [1, 2], 0], #conforming discrete gradient
["L-clamped-Plate" , 0, [1, 2], 0] #nonconforming discrete gradient
]
task_id = 1;
parameterFile = scenarios[task_id][0]
conforming_DiscreteJacobian = scenarios[task_id][1]
gridLevels = scenarios[task_id][2]
INSTRUMENTED = scenarios[task_id][3]
if conforming_DiscreteJacobian:
conformity = "_conforming"
else:
conformity = "_nonconforming"
# x = sys.argv[0]
# print('sys.argv[0]', sys.argv[0])
# print('sys.argv[1]', sys.argv[1])
# print('sys.argv[2]', sys.argv[2])
#
# TAURUS = int(sys.argv[1])
# slurm_array_task_id = int(sys.argv[2])
# TAURUS = 0
# INSTRUMENTED = scenarios[slurm_array_task_id][4]
# print('INSTRUMENTED:', INSTRUMENTED)
#
# print('scenarios[slurm_array_task_id][0]:', scenarios[slurm_array_task_id][0])
# parameterFile = scenarios[slurm_array_task_id][0]
# order = scenarios[slurm_array_task_id][1]
# interpolationMethod = scenarios[slurm_array_task_id][2]
# gridLevels = scenarios[slurm_array_task_id][3]
#Path for parameterFile
# if TAURUS:
# #--- Taurus version
# pythonPath = "/dune/dune-gfe/problems"
# # instrumentedPath = "/dune/dune-gfe/instrumented"
# # resultPath = "/dune/dune-gfe/outputs"
# resultPath = "outputs" + "_" + scenarios[slurm_array_task_id][0]
# instrumentedPath = resultPath + "/instrumented" + "_order" + str(order) + "_" + interpolationMethod
# executablePath = "/dune/dune-gfe/build-cmake/src"
# try:
# os.mkdir(resultPath)
# os.mkdir(instrumentedPath)
# os.mkdir(instrumentedPath + "/mgHistory")
# except OSError as error:
# print(error)
# else :
#--- Local version
pythonPath = "/home/klaus/Desktop/Dune_release/dune-gfe/problems"
# instrumentedPath = '/home/klaus/Desktop/harmonicmapBenchmark/dune-gfe/instrumented'
resultPath = '/home/klaus/Desktop/Dune_release/dune-gfe/outputs' + "_" + scenarios[task_id][0]
instrumentedPath = resultPath + "/instrumented"
executablePath = '/home/klaus/Desktop/Dune_release/dune-gfe/build-cmake/src'
try:
os.mkdir(resultPath)
os.mkdir(instrumentedPath)
os.mkdir(instrumentedPath + "/mgHistory")
except OSError as error:
print(error)
# gridLevels = [1 , 4];
# numReferenceLevels = 6 #grid-level to compute analytical solution
# default = False
# if TAURUS:
# #--- Taurus version
# print('scenarios[slurm_array_task_id][0]:', scenarios[slurm_array_task_id][0])
# parameterFile = scenarios[slurm_array_task_id][0]
# order = scenarios[slurm_array_task_id][1]
# interpolationMethod = scenarios[slurm_array_task_id][2]
# gridLevels = scenarios[slurm_array_task_id][3]
# else :
# #--- Pick experiment manually:
# # EXPERIMENT = 1 #slurm_array_task_id
# print('scenarios[EXPERIMENT][0]:', scenarios[EXPERIMENT][0])
# parameterFile = scenarios[EXPERIMENT][0]
# order = scenarios[EXPERIMENT][1]
# interpolationMethod = scenarios[EXPERIMENT][2]
# gridLevels = scenarios[EXPERIMENT][3]
sys.path.insert(0, pythonPath)
__import__(parameterFile, fromlist=['parameterSet'])
imported_module = importlib.import_module(parameterFile)
parameterSet = getattr(imported_module, "parameterSet")
#--- access parameter settings like:
print('parameterSet.tolerance :', parameterSet.tolerance )
print('parameterSet.structuredGrid :', parameterSet.structuredGrid )
print('parameterSet.numLevels:', parameterSet.numLevels )
dim = parameterSet.dim
# targetDim = parameterSet.targetDim
maxTrustRegionSteps = parameterSet.maxTrustRegionSteps
#############################################
# Compute discrete harmonic maps
#############################################
compute_BendingIsometry = True
# compute_BendingIsometry = False
#############################################
# Compute discretization Error
#############################################
compute_Error = True
# compute_Error = False
# executable = '../build-cmake/src/harmonicmaps-' + str(dim) + "d"
executable = executablePath + "/bending-isometries"
path = os.getcwd()
print("Path: ", path)
############################################
# for order in order_list:
#--- Setup Output-Table:
x = PrettyTable()
constraintError = 0.0
x.title = parameterFile
x.field_names = ["r", "#Triang/#DOF", "# TR-Steps","$E^{Di}(0)$", "$E^{Di}$", "$\delta_1[u_h]$" , "wall-time" ]
x.align["k"] = "l" # Left align city names
x.padding_width = 1 # One space between column edges and contents (default)
rows = []
rows.append(["k", "#Triang/#DOF", "# TR-Steps","$E[u_h^0]$", "$E[u_h^\infty]$", "$\delta_1[u_h]$" , "wall-time" ])
rows_1 = []
rows_2 = []
rows_1.append(["k", "#Triang/#DOF", "$E[u_h^0]$", "$E[u_h^\infty]$"])
rows_2.append(["k", "#Triang/#DOF", "# TR-Steps","$\delta_1[u_h]$" ,"wall-time" ])
if dim == 2:
elements = [8,32,128,512,2048,8192,32768,131072] #Number of elements in Triangulation in 2D
else:
elements = [48,384,3072,24576,196608,1572864] #Number of elements in Triangulation in 3D
# COMPUTE BENDING ISOMETRIES
if compute_BendingIsometry:
computeBendingIsometry(executable,parameterFile,pythonPath,resultPath,instrumentedPath,gridLevels,conforming_DiscreteJacobian ,INSTRUMENTED)
# computeHarmonicMap(minLevel,maxLevel,dim,targetDim,order,interpolationMethod,maxTrustRegionSteps,randomInitialIterate,readConfiguration,configurationFile,perturbation,perturbationRadius,epsilon,tolerance,executable,parameterFile,pythonPath,resultPath)
# --------------------
if compute_Error:
print("-------------------------------------------------------")
print("Now measuring errors with discretizationErrorMode:" + str(parameterSet.computeDiscError.discretizationErrorMode))
# if parameterSet.computeDiscError.discretizationErrorMode=='analytical':
# print("referenceSolution used:", referenceSolution)
# subprocess.call(["echo", "Now measuring errors"])
# EOC_l2_list = ["-","-"]
energyList = []
initial_energyList = []
stepList = []
ndofList = []
timeList = []
# stepTimeList =[]
# EOC_l2_list = []
#
for level in range(gridLevels[0],gridLevels[1]+1):
print("level:", level)
# Measure the discretization errors against the solution on the finest grid
# LOGFILE = "./compute-disc-error_" + str(order) + "_" + str(level) + ".log"
#subprocess.Popen(["../build-cmake/src/compute-disc-error", "compute-disc-error-skyrmions-hexagon.parset",
#"-order", str(order),
#"-level", str(level),
#"-numReferenceLevels", str(maxLevel),
#"-simulationData", "harmonicmaps-result-" + str(order) + "-" + str(level) + ".data",
#"-referenceData", "harmonicmaps-result-" + str(order) + "-" + str(maxLevel) + ".data"])
# ------------ Get Energy / nDofs/Steps etc. -----------------------------------------------------------------------------
LOGFILE_comp = resultPath + "/" + parameterFile + conformity + "_level" + str(level) + ".log"
print('LOGFILE_comp :',LOGFILE_comp )
# LOGFILE_comp = "./harmonicmapsR"+ str(dim) +"_intoR"+ str(targetDim) + "_deg" + str(kappa) + "_"+ interpolationMethod + "_" + str(order) + "_" + str(level) + ".log"
# Read Energy Values:
with open(LOGFILE_comp, 'r') as file:
output = file.read()
try:
# tmp_energy = re.search(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output).group(1) # muss nun nichtmehr am Zeilenanfang stehen! :)
# tmp_energy = re.findall(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output).group(1)
# tmp_energy = re.findall(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
# tmp_energy = re.findall(r'(?m)energy: (-?\d\d?\d?\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
tmp_energy = re.findall(r'(?m)energy: (-?\d?\d?\d?\.?\d+[Ee]?[+\-]?\d?\d?)',output)
tmp_step = re.findall(r'(?m)Step Number: (\d+)',output)
tmp_dofs = re.findall(r'(?m)powerBasis: (\d+)',output)
# tmp_time = re.findall(r'(?m)Time: (-?\d\d?\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
# tmp_time = re.findall(r'(?m)Time: (-?\d+\.\d+[Ee]?[+\-]?\d\d?)',output)
tmp_time = re.findall(r'(?m)Time: (-?\d+\.?\d+\d?)',output)
# tmp_stepTimeList = re.findall(r'(?m)iteration took (-?\d+\.\d+\d?)',output).group(1)
# tmp_stepTimeList = re.findall(r'(?m)Accumulated Timer: (-?\d+\.?\d+\d?)',output)
except AttributeError:
# tmp_energy = re.search(r'(?m)energy: (-?\d\d?\.\d+[Ee]?[+\-]?\d\d?)',output) # muss nun nichtmehr am Zeilenanfang stehen! :)
tmp_energy = re.search(r'(?m)energy: (-?\d?\d?\d?\.?\d+[Ee]?[+\-]?\d?\d?)',output)
tmp_step = re.findall(r'(?m)Step Number: (\d+)',output)
tmp_dofs = re.findall(r'(?m)powerBasis: (\d+)',output)
# tmp_time = re.findall(r'(?m)Time: (-?\d\d?\d?\d?\d?\.\d+[Ee]?[+\-]?\d\d?)',output)
# tmp_time = re.findall(r'(?m)Time: (-?\d+\.\d+[Ee]?[+\-]?\d\d?)',output)
tmp_time = re.findall(r'(?m)Time: (-?\d+\.?\d+\d?)',output)
# tmp_stepTimeList = re.findall(r'(?m)iteration took (-?\d+\.\d+\d?)',output)
# tmp_stepTimeList = re.findall(r'(?m)Accumulated Timer: (-?\d+\.?\d+\d?)',output)
if maxTrustRegionSteps>0 :
print('tmp_step:',tmp_step)
print('tmp_step last:', int(tmp_step[-1]))
# print('tmp_dofs:',tmp_dofs)
# print('tmp_dofs last:', int(tmp_dofs[-1]))
print('type tmp_energy:', type(tmp_energy))
print('tmp_energy:', tmp_energy)
print('tmp_energy last:', tmp_energy[-1])
print('tmp_time:', tmp_time)
# print('tmp_stepTimeList', tmp_stepTimeList)
# stepTimeList = [float(i) for i in tmp_stepTimeList] # Convert list to float
# print('stepTimeList', stepTimeList)
stepEnergyList = [float(i) for i in tmp_energy] # Convert list to float
print('stepEnergyList', stepEnergyList)
energy = float(tmp_energy[-1]) #[1] since otherwise it recognizes "2" from L2error...
energyList.append(energy)
initialEnergy = float(tmp_energy[0])
initial_energyList.append(initialEnergy)
steps = int(tmp_step[-1])+1 #starts with zero therefore +1
stepList.append(steps)
# ndofs = int(tmp_dofs[-1])
# ndofList.append(ndofs)
time = float(tmp_time[-1])
timeList.append(time)
#######################################################################
print('--------------------- GridLevel:' + str(level) + '----------------------')
# print('(accumulated) wall-time for each iteration-step', stepTimeList)
print('energy for each iteration step', stepEnergyList)
print('----------------------------------------------------------------------------')
# -----------------------------------------------------------------------------------------------
EOC_l2 = "-"
EOC_h1 = "-"
EOC_energy = "-"
L2_fine = "-"
L2_coarse = "-"
H1_fine = "-"
H1_coarse = "-"
# constraintError = "-"
if maxTrustRegionSteps > 0 :
print('energyList', energyList)
# currentEnergy = energyList[level-1]
currentEnergy = energyList[-1] # better like this?
print('currentEnergy', currentEnergy)
# currentDofs = ndofList[level-1]
# currentDofs = ndofList[-1]
# print('currentDofs', currentDofs)
# currentSteps = stepList[level-1]
currentSteps = stepList[-1]
print('currentSteps', currentSteps)
print('initial_energyList', initial_energyList)
currentInitialEnergy = initial_energyList[-1]
print('currentInitialEnergy', currentInitialEnergy)
else:
currentEnergy = "-"
currentDofs = "-"
currentSteps = "-"
currentInitialEnergy = "-"
print('timeList:', timeList)
currentTime = timeList[-1]
print('currentSteps:', currentSteps)
x.add_row([level, elements[level-1], currentSteps, currentInitialEnergy, currentEnergy,constraintError, currentTime])
# rows.append([level, currentDofs, str(currentSteps), currentInitialEnergy, currentEnergy,constraintError, EOC_l2, EOC_h1, EOC_energy , currentTime])
rows.append([level, elements[level-1], currentSteps,r"\num{" +str(currentInitialEnergy)+"}", "\\"+"num{"+str(currentEnergy)+"}","\\num{"+str(constraintError)+"}" , currentTime]) #new: add \num{ }
rows_1.append([level, elements[level-1], r"\num{" +str(currentInitialEnergy)+"}", "\\"+"num{"+str(currentEnergy)+"}"]) #new: add \num{ }
rows_2.append([level, elements[level-1], currentSteps,"\\num{"+str(constraintError)+"}", currentTime]) #new: add \num{ }
## Add extra column:
# print(*EOC_l2_list)
# x.add_column('EOC', EOC_l2_list)
# Write Table to Text-File:
# if pseudoRandomInitialIterate:
# tablefile = open("EOC-table_R"+str(dim)+"intoR"+ str(targetDim)+ "_PRNG" + "_" + interpolationMethod+ "_order" +str(order) + "tol" + str(tolerance) + ".txt", "w")
# else:
# tablefile = open("EOC-table_R"+str(dim)+"intoR"+ str(targetDim)+ "_deg" + str(kappa) + "_" + interpolationMethod+ "_order" +str(order) + "tol" + str(tolerance) + ".txt", "w")
# if pseudoRandomInitialIterate:
# tablefile = open(parameterFile + "_order" + str(order) + "_" + interpolationMethod + "_PRNG" + "_table" + ".txt", "w")
# else:
tablefile = open(resultPath + "/" + parameterFile + conformity + "_table" + ".txt", "w")
tablefile.write(str(x))
tablefile.write('\n')
tablefile.write(str(tabulate(rows_1, headers='firstrow', tablefmt='latex_raw')))
tablefile.write(str(tabulate(rows_2, headers='firstrow', tablefmt='latex_raw')))
tablefile.close()
# print (pretty)Table
print(x)
print(tabulate(rows_1, headers='firstrow', tablefmt='latex_raw'))
print(tabulate(rows_2, headers='firstrow', tablefmt='latex_raw'))
########## end of kappa-loop ##########################