From fa92caaef90e64fd858bde336015dea0a2a63dcc Mon Sep 17 00:00:00 2001
From: Lars Bilke <lars.bilke@ufz.de>
Date: Thu, 15 Oct 2020 09:06:21 +0200
Subject: [PATCH] Format python: Tests/.
---
.../cube_1x1x1_SteadyStateDiffusion/cube.py | 122 ++++---
.../bcs_laplace_eq.py | 15 +-
.../sin_x_sin_y_source_term.py | 19 +-
.../square_1e1_neumann-insitu.py | 80 +++--
...mn_nonequilibrium_sigma_generate_values.py | 14 +-
.../PythonHertzContact/gen-unit-circle.py | 39 ++-
.../PythonHertzContact/hertz_contact_bc.py | 46 +--
.../Linear/PythonHertzContact/post.py | 127 ++++---
.../Mechanics/Linear/PythonPiston/chamber.py | 14 +-
.../Linear/PythonPiston/piston_bc.py | 5 +-
.../Mechanics/Linear/PythonPiston/post.py | 21 +-
Tests/Data/Parabolic/T/1D_neumann/plotLine.py | 325 +++++++++---------
.../T/1D_neumann/temperature_analytical.py | 12 +-
.../Parabolic/T/3D_3BHEs_array/bcs_tespy.py | 105 +++---
.../T/3D_3BHEs_array/bcs_tespy_closedloop.py | 105 +++---
.../Parabolic/T/3D_3BHEs_array/pre/3bhes.py | 168 ++++++---
.../T/3D_3BHEs_array/pre/3bhes_closedloop.py | 151 +++++---
.../Data/ThermoMechanics/BDT/generate_ref.py | 6 +-
18 files changed, 806 insertions(+), 568 deletions(-)
diff --git a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube.py b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube.py
index 87325b914ff..cc00ae7edc0 100644
--- a/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube.py
+++ b/Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube.py
@@ -1,76 +1,96 @@
-
from paraview.simple import *
from paraview import coprocessing
-#--------------------------------------------------------------
+# --------------------------------------------------------------
# Code generated from cpstate.py to create the CoProcessor.
# ParaView 5.2.0 64 bits
# ----------------------- CoProcessor definition -----------------------
+
def CreateCoProcessor():
- def _CreatePipeline(coprocessor, datadescription):
- class Pipeline:
- # state file generated using paraview version 5.2.0
-
- # ----------------------------------------------------------------
- # setup the data processing pipelines
- # ----------------------------------------------------------------
-
- #### disable automatic camera reset on 'Show'
- paraview.simple._DisableFirstRenderCameraReset()
-
- # create a new 'XML Unstructured Grid Reader'
- # create a producer from a simulation input
- cube_1e0_pcs_0_ts_0_t_0000000vtu = coprocessor.CreateProducer(datadescription, 'input')
-
- # create a new 'Contour'
- contour1 = Contour(Input=cube_1e0_pcs_0_ts_0_t_0000000vtu)
- contour1.ContourBy = ['POINTS', 'D1_left_front_N1_right']
- contour1.ComputeScalars = 1
- contour1.Isosurfaces = [1.0, 1.0750344444444444, 1.1500688888888888, 1.2251033333333334, 1.3001377777777778, 1.3751722222222222, 1.4502066666666669, 1.5252411111111113, 1.6002755555555557, 1.67531]
- contour1.PointMergeMethod = 'Uniform Binning'
-
- # create a new 'Parallel PolyData Writer'
- parallelPolyDataWriter1 = servermanager.writers.XMLPPolyDataWriter(Input=contour1)
-
- # register the writer with coprocessor
- # and provide it with information such as the filename to use,
- # how frequently to write the data, etc.
- coprocessor.RegisterWriter(parallelPolyDataWriter1, filename='cube_1e1_%t.pvtp', freq=1)
-
- # ----------------------------------------------------------------
- # finally, restore active source
- SetActiveSource(contour1)
- # ----------------------------------------------------------------
- return Pipeline()
-
- class CoProcessor(coprocessing.CoProcessor):
- def CreatePipeline(self, datadescription):
- self.Pipeline = _CreatePipeline(self, datadescription)
-
- coprocessor = CoProcessor()
- # these are the frequencies at which the coprocessor updates.
- freqs = {'input': [1]}
- coprocessor.SetUpdateFrequencies(freqs)
- return coprocessor
-
-#--------------------------------------------------------------
+ def _CreatePipeline(coprocessor, datadescription):
+ class Pipeline:
+ # state file generated using paraview version 5.2.0
+
+ # ----------------------------------------------------------------
+ # setup the data processing pipelines
+ # ----------------------------------------------------------------
+
+ #### disable automatic camera reset on 'Show'
+ paraview.simple._DisableFirstRenderCameraReset()
+
+ # create a new 'XML Unstructured Grid Reader'
+ # create a producer from a simulation input
+ cube_1e0_pcs_0_ts_0_t_0000000vtu = coprocessor.CreateProducer(
+ datadescription, "input"
+ )
+
+ # create a new 'Contour'
+ contour1 = Contour(Input=cube_1e0_pcs_0_ts_0_t_0000000vtu)
+ contour1.ContourBy = ["POINTS", "D1_left_front_N1_right"]
+ contour1.ComputeScalars = 1
+ contour1.Isosurfaces = [
+ 1.0,
+ 1.0750344444444444,
+ 1.1500688888888888,
+ 1.2251033333333334,
+ 1.3001377777777778,
+ 1.3751722222222222,
+ 1.4502066666666669,
+ 1.5252411111111113,
+ 1.6002755555555557,
+ 1.67531,
+ ]
+ contour1.PointMergeMethod = "Uniform Binning"
+
+ # create a new 'Parallel PolyData Writer'
+ parallelPolyDataWriter1 = servermanager.writers.XMLPPolyDataWriter(
+ Input=contour1
+ )
+
+ # register the writer with coprocessor
+ # and provide it with information such as the filename to use,
+ # how frequently to write the data, etc.
+ coprocessor.RegisterWriter(
+ parallelPolyDataWriter1, filename="cube_1e1_%t.pvtp", freq=1
+ )
+
+ # ----------------------------------------------------------------
+ # finally, restore active source
+ SetActiveSource(contour1)
+ # ----------------------------------------------------------------
+
+ return Pipeline()
+
+ class CoProcessor(coprocessing.CoProcessor):
+ def CreatePipeline(self, datadescription):
+ self.Pipeline = _CreatePipeline(self, datadescription)
+
+ coprocessor = CoProcessor()
+ # these are the frequencies at which the coprocessor updates.
+ freqs = {"input": [1]}
+ coprocessor.SetUpdateFrequencies(freqs)
+ return coprocessor
+
+
+# --------------------------------------------------------------
# Global variables that will hold the pipeline for each timestep
# Creating the CoProcessor object, doesn't actually create the ParaView pipeline.
# It will be automatically setup when coprocessor.UpdateProducers() is called the
# first time.
coprocessor = CreateCoProcessor()
-#--------------------------------------------------------------
+# --------------------------------------------------------------
# Enable Live-Visualizaton with ParaView
coprocessor.EnableLiveVisualization(False, 1)
# ---------------------- Data Selection method ----------------------
+
def RequestDataDescription(datadescription):
"Callback to populate the request for current timestep"
global coprocessor
@@ -86,8 +106,10 @@ def RequestDataDescription(datadescription):
# pipeline.
coprocessor.LoadRequestedData(datadescription)
+
# ------------------------ Processing method ------------------------
+
def DoCoProcessing(datadescription):
"Callback to do co-processing for current timestep"
global coprocessor
diff --git a/Tests/Data/Elliptic/square_1x1_SteadyStateDiffusion_Python/bcs_laplace_eq.py b/Tests/Data/Elliptic/square_1x1_SteadyStateDiffusion_Python/bcs_laplace_eq.py
index e6d5c8dd7d8..89c5c4e8e5c 100644
--- a/Tests/Data/Elliptic/square_1x1_SteadyStateDiffusion_Python/bcs_laplace_eq.py
+++ b/Tests/Data/Elliptic/square_1x1_SteadyStateDiffusion_Python/bcs_laplace_eq.py
@@ -1,17 +1,17 @@
-
import OpenGeoSys
from math import pi, sin, cos, sinh, cosh
-a = 2.0*pi/3.0
+a = 2.0 * pi / 3.0
# analytical solution used to set the Dirichlet BCs
def solution(x, y):
- return sin(a*x) * sinh(a*y)
+ return sin(a * x) * sinh(a * y)
+
# gradient of the analytical solution used to set the Neumann BCs
def grad_solution(x, y):
- return a * cos(a*x) * sinh(a*y), \
- a * sin(a*x) * cosh(a*y)
+ return a * cos(a * x) * sinh(a * y), a * sin(a * x) * cosh(a * y)
+
# Dirichlet BCs
class BCTop(OpenGeoSys.BoundaryCondition):
@@ -21,6 +21,7 @@ class BCTop(OpenGeoSys.BoundaryCondition):
value = solution(x, y)
return (True, value)
+
class BCLeft(OpenGeoSys.BoundaryCondition):
def getDirichletBCValue(self, t, coords, node_id, primary_vars):
x, y, z = coords
@@ -28,6 +29,7 @@ class BCLeft(OpenGeoSys.BoundaryCondition):
value = solution(x, y)
return (True, value)
+
class BCBottom(OpenGeoSys.BoundaryCondition):
def getDirichletBCValue(self, t, coords, node_id, primary_vars):
x, y, z = coords
@@ -35,13 +37,14 @@ class BCBottom(OpenGeoSys.BoundaryCondition):
value = solution(x, y)
return (True, value)
+
# Neumann BC
class BCRight(OpenGeoSys.BoundaryCondition):
def getFlux(self, t, coords, primary_vars):
x, y, z = coords
assert x == 1.0 and z == 0.0
value = grad_solution(x, y)[0]
- Jac = [ 0.0 ] # value does not depend on primary variable
+ Jac = [0.0] # value does not depend on primary variable
return (True, value, Jac)
diff --git a/Tests/Data/Elliptic/square_1x1_SteadyStateDiffusion_Python/sin_x_sin_y_source_term.py b/Tests/Data/Elliptic/square_1x1_SteadyStateDiffusion_Python/sin_x_sin_y_source_term.py
index 60766042860..310bb0ff18e 100644
--- a/Tests/Data/Elliptic/square_1x1_SteadyStateDiffusion_Python/sin_x_sin_y_source_term.py
+++ b/Tests/Data/Elliptic/square_1x1_SteadyStateDiffusion_Python/sin_x_sin_y_source_term.py
@@ -1,24 +1,29 @@
-
import OpenGeoSys
from math import pi, sin
-a = 2.0*pi
-b = 2.0*pi
+a = 2.0 * pi
+b = 2.0 * pi
+
def solution(x, y):
- return sin(a*x-pi/2.0) * sin(b*y-pi/2.0)
+ return sin(a * x - pi / 2.0) * sin(b * y - pi / 2.0)
+
# - laplace(solution) = source term
def laplace_solution(x, y):
- return a*a * sin(a*x-pi/2.0) * sin(b*y-pi/2.0) + b*b * sin(a*x-pi/2.0) * sin(b*y-pi/2.0)
+ return a * a * sin(a * x - pi / 2.0) * sin(b * y - pi / 2.0) + b * b * sin(
+ a * x - pi / 2.0
+ ) * sin(b * y - pi / 2.0)
+
# source term for the benchmark
class SinXSinYSourceTerm(OpenGeoSys.SourceTerm):
def getFlux(self, t, coords, primary_vars):
x, y, z = coords
- value = laplace_solution(x,y)
- Jac = [ 0.0 ]
+ value = laplace_solution(x, y)
+ Jac = [0.0]
return (value, Jac)
+
# instantiate source term object referenced in OpenGeoSys' prj file
sinx_siny_source_term = SinXSinYSourceTerm()
diff --git a/Tests/Data/EllipticPETSc/square_1e1_neumann-insitu.py b/Tests/Data/EllipticPETSc/square_1e1_neumann-insitu.py
index a1073072bd5..7fb0cee4588 100644
--- a/Tests/Data/EllipticPETSc/square_1e1_neumann-insitu.py
+++ b/Tests/Data/EllipticPETSc/square_1e1_neumann-insitu.py
@@ -1,4 +1,4 @@
-#--------------------------------------------------------------
+# --------------------------------------------------------------
# Global timestep output options
timeStepToStartOutputAt = 0
@@ -12,15 +12,15 @@ rescale_lookuptable = False
requestSpecificArrays = False
# a root directory under which all Catalyst output goes
-rootDirectory = ''
+rootDirectory = ""
# makes a cinema D index table
make_cinema_table = False
-#--------------------------------------------------------------
+# --------------------------------------------------------------
# Code generated from cpstate.py to create the CoProcessor.
# paraview version 5.8.0
-#--------------------------------------------------------------
+# --------------------------------------------------------------
from paraview.simple import *
from paraview import coprocessing
@@ -46,13 +46,13 @@ def CreateCoProcessor():
paraview.simple._DisableFirstRenderCameraReset()
# create a new 'Superquadric'
- ogs_output = coprocessor.CreateProducer(datadescription, 'input')
+ ogs_output = coprocessor.CreateProducer(datadescription, "input")
# create a new 'Clip'
clip1 = Clip(Input=ogs_output)
- clip1.ClipType = 'Plane'
- clip1.HyperTreeGridClipper = 'Plane'
- clip1.Scalars = ['POINTS', 'D1_left_bottom_N1_right']
+ clip1.ClipType = "Plane"
+ clip1.HyperTreeGridClipper = "Plane"
+ clip1.Scalars = ["POINTS", "D1_left_bottom_N1_right"]
clip1.Value = 1.3376572416618493
# init the 'Plane' selected for 'ClipType'
@@ -63,13 +63,20 @@ def CreateCoProcessor():
# create a new 'Contour'
contour1 = Contour(Input=ogs_output)
- contour1.ContourBy = ['POINTS', 'D1_left_bottom_N1_right']
+ contour1.ContourBy = ["POINTS", "D1_left_bottom_N1_right"]
contour1.Isosurfaces = [
- 1.0, 1.075034942591522, 1.1500698851830442, 1.2251048277745662,
- 1.3001397703660882, 1.3751747129576102, 1.4502096555491324,
- 1.5252445981406544, 1.6002795407321764, 1.6753144833236986
+ 1.0,
+ 1.075034942591522,
+ 1.1500698851830442,
+ 1.2251048277745662,
+ 1.3001397703660882,
+ 1.3751747129576102,
+ 1.4502096555491324,
+ 1.5252445981406544,
+ 1.6002795407321764,
+ 1.6753144833236986,
]
- contour1.PointMergeMethod = 'Uniform Binning'
+ contour1.PointMergeMethod = "Uniform Binning"
# ----------------------------------------------------------------
# finally, restore active source
@@ -77,31 +84,35 @@ def CreateCoProcessor():
# ----------------------------------------------------------------
# Now any catalyst writers
- xMLPUnstructuredGridWriter1 = servermanager.writers.XMLPUnstructuredGridWriter(
- Input=clip1)
+ xMLPUnstructuredGridWriter1 = (
+ servermanager.writers.XMLPUnstructuredGridWriter(Input=clip1)
+ )
coprocessor.RegisterWriter(
xMLPUnstructuredGridWriter1,
- filename='square_1e1_neumann_clip1_%t.pvtu',
+ filename="square_1e1_neumann_clip1_%t.pvtu",
freq=1,
paddingamount=0,
- DataMode='Appended',
- HeaderType='UInt64',
+ DataMode="Appended",
+ HeaderType="UInt64",
EncodeAppendedData=False,
- CompressorType='None',
- CompressionLevel='6')
+ CompressorType="None",
+ CompressionLevel="6",
+ )
xMLPPolyDataWriter1 = servermanager.writers.XMLPPolyDataWriter(
- Input=contour1)
+ Input=contour1
+ )
coprocessor.RegisterWriter(
xMLPPolyDataWriter1,
- filename='square_1e1_neumann_contour1_%t.pvtp',
+ filename="square_1e1_neumann_contour1_%t.pvtp",
freq=1,
paddingamount=0,
- DataMode='Appended',
- HeaderType='UInt64',
+ DataMode="Appended",
+ HeaderType="UInt64",
EncodeAppendedData=False,
- CompressorType='None',
- CompressionLevel='6')
+ CompressorType="None",
+ CompressionLevel="6",
+ )
return Pipeline()
@@ -113,8 +124,7 @@ def CreateCoProcessor():
# these are the frequencies at which the coprocessor updates.
freqs = {}
coprocessor.SetUpdateFrequencies(freqs)
- coprocessor.SetInitialOutputOptions(timeStepToStartOutputAt,
- forceOutputAtFirstCall)
+ coprocessor.SetInitialOutputOptions(timeStepToStartOutputAt, forceOutputAtFirstCall)
if rootDirectory:
coprocessor.SetRootDirectory(rootDirectory)
@@ -125,14 +135,14 @@ def CreateCoProcessor():
return coprocessor
-#--------------------------------------------------------------
+# --------------------------------------------------------------
# Global variable that will hold the pipeline for each timestep
# Creating the CoProcessor object, doesn't actually create the ParaView pipeline.
# It will be automatically setup when coprocessor.UpdateProducers() is called the
# first time.
coprocessor = CreateCoProcessor()
-#--------------------------------------------------------------
+# --------------------------------------------------------------
# Enable Live-Visualizaton with ParaView and the update frequency
coprocessor.EnableLiveVisualization(False, 1)
@@ -163,10 +173,12 @@ def DoCoProcessing(datadescription):
coprocessor.WriteData(datadescription)
# Write image capture (Last arg: rescale lookup table), if appropriate.
- coprocessor.WriteImages(datadescription,
- rescale_lookuptable=rescale_lookuptable,
- image_quality=0,
- padding_amount=imageFileNamePadding)
+ coprocessor.WriteImages(
+ datadescription,
+ rescale_lookuptable=rescale_lookuptable,
+ image_quality=0,
+ padding_amount=imageFileNamePadding,
+ )
# Live Visualization, if enabled.
coprocessor.DoLiveVisualization(datadescription, "localhost", 22222)
diff --git a/Tests/Data/Mechanics/InitialStates/soil_column_nonequilibrium_sigma_generate_values.py b/Tests/Data/Mechanics/InitialStates/soil_column_nonequilibrium_sigma_generate_values.py
index db5b66ac85f..1144470e77d 100755
--- a/Tests/Data/Mechanics/InitialStates/soil_column_nonequilibrium_sigma_generate_values.py
+++ b/Tests/Data/Mechanics/InitialStates/soil_column_nonequilibrium_sigma_generate_values.py
@@ -8,10 +8,11 @@ def createStressArray(points):
sigma.SetNumberOfTuples(points.GetNumberOfPoints())
for i in range(points.GetNumberOfPoints()):
y = points.GetPoint(i)[1]
- sigma.SetTuple4(i, 2200 * 9.81 * 0.8 * (y - 100),
- -0.5 * 1e6 + 2200 * 9.81 * (y - 100), 0, 0)
- #sigma.SetTuple4(i, 4*i, 4*i + 1, 4*i + 2, 4*i + 3) # for debugging
- sigma.SetName('nonequilibrium_stress')
+ sigma.SetTuple4(
+ i, 2200 * 9.81 * 0.8 * (y - 100), -0.5 * 1e6 + 2200 * 9.81 * (y - 100), 0, 0
+ )
+ # sigma.SetTuple4(i, 4*i, 4*i + 1, 4*i + 2, 4*i + 3) # for debugging
+ sigma.SetName("nonequilibrium_stress")
return sigma
@@ -32,14 +33,15 @@ def writeDataToFile(mesh, filename):
w.SetInputData(mesh)
w.Write()
+
def addStressArraysToMesh(mesh):
mesh.GetPointData().AddArray(createStressArrayForNodes(mesh))
mesh.GetCellData().AddArray(createStressArrayForCells(mesh))
r = vtkXMLUnstructuredGridReader()
-r.SetFileName('soil_column.vtu')
+r.SetFileName("soil_column.vtu")
r.Update()
m = r.GetOutput()
addStressArraysToMesh(m)
-writeDataToFile(m, 'soil_column_nonequilibrium_sigma.vtu')
+writeDataToFile(m, "soil_column_nonequilibrium_sigma.vtu")
diff --git a/Tests/Data/Mechanics/Linear/PythonHertzContact/gen-unit-circle.py b/Tests/Data/Mechanics/Linear/PythonHertzContact/gen-unit-circle.py
index ef4bf9072ae..a4a20cdfbcc 100755
--- a/Tests/Data/Mechanics/Linear/PythonHertzContact/gen-unit-circle.py
+++ b/Tests/Data/Mechanics/Linear/PythonHertzContact/gen-unit-circle.py
@@ -12,10 +12,10 @@ def distribute_points_evenly(c2):
assert np.sqrt(c2.shape[0]).is_integer()
CELLS_PER_DIRECTION = int(np.sqrt(c2.shape[0])) - 1
- r2 = np.sqrt(c2[:,0]**2 + c2[:,1]**2)
- alpha2 = np.arctan2(c2[:,1], c2[:,0])
+ r2 = np.sqrt(c2[:, 0] ** 2 + c2[:, 1] ** 2)
+ alpha2 = np.arctan2(c2[:, 1], c2[:, 0])
- bins = [ [] for _ in range(CELLS_PER_DIRECTION+1) ]
+ bins = [[] for _ in range(CELLS_PER_DIRECTION + 1)]
nbin = np.round(r2 * CELLS_PER_DIRECTION).astype(int)
for node, b in enumerate(nbin):
bins[b].append(node)
@@ -53,10 +53,10 @@ a = np.empty(coords.shape[0])
for i, (x, y, z) in enumerate(coords):
if x > y:
- R = np.sqrt(1 + (y/x)**2)
+ R = np.sqrt(1 + (y / x) ** 2)
a[i] = 1.0 / R
elif x < y:
- R = np.sqrt(1 + (x/y)**2)
+ R = np.sqrt(1 + (x / y) ** 2)
a[i] = 1.0 / R
else:
a[i] = np.sqrt(0.5)
@@ -81,28 +81,32 @@ writer.SetInputData(grid)
writer.Write()
# extract boundary of unit circle
-R_squared = new_coords[:,0]**2 + new_coords[:,1]**2
-phi = np.arctan2(new_coords[:,1], new_coords[:,0])
+R_squared = new_coords[:, 0] ** 2 + new_coords[:, 1] ** 2
+phi = np.arctan2(new_coords[:, 1], new_coords[:, 0])
idcs = np.where(abs(R_squared - 1) < 1e-8)[0]
idcs = idcs[np.argsort(phi[idcs])] # sorted with ascending polar angle
with open(out_geom, "w") as fh:
- fh.write("""<?xml version="1.0" encoding="ISO-8859-1"?>
+ fh.write(
+ """<?xml version="1.0" encoding="ISO-8859-1"?>
<OpenGeoSysGLI>
<name>geom</name>
<points>
<point id="0" x="0" y="0" z="0" name="center"/>
<point id="1" x="0" y="1" z="0" name="top"/>
<point id="2" x="1" y="0" z="0"/>
-""")
+"""
+ )
for i, (x, y, z) in enumerate(new_coords[idcs]):
- fh.write(' <point id="{}" x="{}" y="{}" z="{}" />\n'.format(
- i+3, x, y, z))
+ fh.write(
+ ' <point id="{}" x="{}" y="{}" z="{}" />\n'.format(i + 3, x, y, z)
+ )
- fh.write("""
+ fh.write(
+ """
</points>
<polylines>
@@ -114,11 +118,14 @@ with open(out_geom, "w") as fh:
<pnt>0</pnt>
<pnt>2</pnt>
</polyline>
- <polyline id="2" name="outer">\n""")
+ <polyline id="2" name="outer">\n"""
+ )
for i in range(len(idcs)):
- fh.write(" <pnt>{}</pnt>\n".format(i+3))
+ fh.write(" <pnt>{}</pnt>\n".format(i + 3))
- fh.write("""</polyline>
+ fh.write(
+ """</polyline>
</polylines>
-</OpenGeoSysGLI>""")
+</OpenGeoSysGLI>"""
+ )
diff --git a/Tests/Data/Mechanics/Linear/PythonHertzContact/hertz_contact_bc.py b/Tests/Data/Mechanics/Linear/PythonHertzContact/hertz_contact_bc.py
index 4abfa50b91b..ba293a82b7b 100644
--- a/Tests/Data/Mechanics/Linear/PythonHertzContact/hertz_contact_bc.py
+++ b/Tests/Data/Mechanics/Linear/PythonHertzContact/hertz_contact_bc.py
@@ -11,26 +11,26 @@ class HertzContactBC(OpenGeoSys.BoundaryCondition):
def __init__(self):
super(HertzContactBC, self).__init__()
- self._first_node = None # ID of the first node of this BC's geometry
+ self._first_node = None # ID of the first node of this BC's geometry
self._t_old = START_TIME - 1.0 # time of previous invocation of this BC
- self._boundary_x_coords = [] # the x coordinates of all boundary nodes
-
+ self._boundary_x_coords = [] # the x coordinates of all boundary nodes
def _init_timestep(self):
"""Initializes the internal state at the beginning of a new timestep."""
- self._a_range = [ 0.0, SPHERE_RADIUS ] # range of possible contact radii
- self._a_est = SPHERE_RADIUS # estimated contact radius
-
- self._max_x_with_y_excess = -1.0 # maximum x value where a node is above the contact line
+ self._a_range = [0.0, SPHERE_RADIUS] # range of possible contact radii
+ self._a_est = SPHERE_RADIUS # estimated contact radius
- self._tendency_is_up = True # whether the contact area is supposed to grow in the current iteration
- self._a_curr = 0.0 # the radius of the contact area in this iteration (continuously updated)
- self._a_prev = 0.0 # the radius of the contact area that was prescripted in the previous iteration
+ self._max_x_with_y_excess = (
+ -1.0
+ ) # maximum x value where a node is above the contact line
- self._iteration = 0 # the nonlinear solver iteration number
+ self._tendency_is_up = True # whether the contact area is supposed to grow in the current iteration
+ self._a_curr = 0.0 # the radius of the contact area in this iteration (continuously updated)
+ self._a_prev = 0.0 # the radius of the contact area that was prescripted in the previous iteration
+ self._iteration = 0 # the nonlinear solver iteration number
def _init_iteration(self):
"""Initializes the internal state at the beginning of a new nonlinear solver iteration."""
@@ -66,9 +66,11 @@ class HertzContactBC(OpenGeoSys.BoundaryCondition):
self._iteration += 1
- print("BC: a_est={:.4f}, a_prev={:.4f} ({:.4f}, {:.4f})".format(
- self._a_est, self._a_prev, self._a_range[0], self._a_range[1]))
-
+ print(
+ "BC: a_est={:.4f}, a_prev={:.4f} ({:.4f}, {:.4f})".format(
+ self._a_est, self._a_prev, self._a_range[0], self._a_range[1]
+ )
+ )
def getDirichletBCValue(self, t, coords, node_id, primary_vars):
if self._t_old < t:
@@ -93,10 +95,15 @@ class HertzContactBC(OpenGeoSys.BoundaryCondition):
try:
# check that we are at the outer boundary
- assert abs(x**2 + y**2 + z**2 - SPHERE_RADIUS**2) < 1e-15
+ assert abs(x ** 2 + y ** 2 + z ** 2 - SPHERE_RADIUS ** 2) < 1e-15
except:
- print("assert abs(x**2 + y**2 + z**2 - 1.0) < 1e-15",
- x, y, z, x**2 + y**2 + z**2 - SPHERE_RADIUS**2)
+ print(
+ "assert abs(x**2 + y**2 + z**2 - 1.0) < 1e-15",
+ x,
+ y,
+ z,
+ x ** 2 + y ** 2 + z ** 2 - SPHERE_RADIUS ** 2,
+ )
raise
y_deformed = y + uy
@@ -128,20 +135,17 @@ class HertzContactBC(OpenGeoSys.BoundaryCondition):
res = (True, y_top - y)
elif self._boundary_x_coords_are_initialized(t):
idx = self._boundary_x_coords.index(x)
- if idx != 0 and self._boundary_x_coords[idx-1] == self._a_prev:
+ if idx != 0 and self._boundary_x_coords[idx - 1] == self._a_prev:
res = (True, y_top - y)
-
return res
-
@staticmethod
def _get_y_top(t):
"""Returns the y-position of the contact area depending on the load step t."""
return 1.0 - 0.005 * t
-
def _boundary_x_coords_are_initialized(self, t):
return self._iteration >= 2
diff --git a/Tests/Data/Mechanics/Linear/PythonHertzContact/post.py b/Tests/Data/Mechanics/Linear/PythonHertzContact/post.py
index 8dc81bf1acb..1e4ffcd2d79 100755
--- a/Tests/Data/Mechanics/Linear/PythonHertzContact/post.py
+++ b/Tests/Data/Mechanics/Linear/PythonHertzContact/post.py
@@ -20,23 +20,22 @@ lambda_ = E_12 * nu_12 / (1 + nu_12) / (1 - 2 * nu_12)
mu = E_12 / 2.0 / (1 + nu_12)
G_12 = mu
-kappa = 0.5 * G_12 / R / (1-nu_12)
+kappa = 0.5 * G_12 / R / (1 - nu_12)
print("kappa:", kappa)
-C = lambda_ * np.matrix([
- 1, 1, 1, 0,
- 1, 1, 1, 0,
- 1, 1, 1, 0,
- 0, 0, 0, 0 ]).reshape(4, 4) \
- + 2 * mu * np.identity(4)
+C = lambda_ * np.matrix([1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0]).reshape(
+ 4, 4
+) + 2 * mu * np.identity(4)
+
def p_contact(r, r_contact):
- return kappa * np.sqrt(r_contact**2 - r**2)
+ return kappa * np.sqrt(r_contact ** 2 - r ** 2)
### helpers ##############################################
import os
+
try:
import xml.etree.cElementTree as ET
except:
@@ -48,6 +47,7 @@ def relpathfrom(origin, relpath):
return relpath
return os.path.join(origin, relpath)
+
def read_pvd_file(fn):
try:
path = fn.name
@@ -56,28 +56,34 @@ def read_pvd_file(fn):
pathroot = os.path.dirname(path)
pvdtree = ET.parse(fn)
node = pvdtree.getroot()
- if node.tag != "VTKFile": return None, None
+ if node.tag != "VTKFile":
+ return None, None
children = list(node)
- if len(children) != 1: return None, None
+ if len(children) != 1:
+ return None, None
node = children[0]
- if node.tag != "Collection": return None, None
+ if node.tag != "Collection":
+ return None, None
ts = []
fs = []
for child in node:
- if child.tag != "DataSet": return None, None
+ if child.tag != "DataSet":
+ return None, None
ts.append(float(child.get("timestep")))
fs.append(relpathfrom(pathroot, child.get("file")))
return ts, fs
+
### helpers end ##########################################
def get_y_top(t):
return 1.0 - 0.005 * t
+
ts, fns = read_pvd_file(pvd_file)
reader = vtk.vtkXMLUnstructuredGridReader()
@@ -131,7 +137,7 @@ for t, fn in zip(ts, fns):
disp_2d = vtk_to_numpy(grid.GetPointData().GetArray("displacement"))
disp_3d = np.zeros((disp_2d.shape[0], 3))
- disp_3d[:,(0,1)] = disp_2d
+ disp_3d[:, (0, 1)] = disp_2d
disp_3d_vtk = numpy_to_vtk(disp_3d, 1)
disp_3d_vtk.SetName("u")
@@ -141,34 +147,35 @@ for t, fn in zip(ts, fns):
if False:
# compute strain
def strain_triangle_axi(cell, point_data, strain_data):
- cell_pts = np.matrix(vtk_to_numpy(cell.GetPoints().GetData())[:,:-1])
+ cell_pts = np.matrix(vtk_to_numpy(cell.GetPoints().GetData())[:, :-1])
assert cell_pts.shape[0] == 3 # triangles
assert point_data.shape[1] == 2 # 2D vector field
- node_ids = [ cell.GetPointId(i) for i in range(cell.GetNumberOfPoints()) ]
+ node_ids = [cell.GetPointId(i) for i in range(cell.GetNumberOfPoints())]
# interpolation using barycentric coordinates on linear triangles
- T = np.matrix(np.empty((2,2)))
- T[:,0] = (cell_pts[0,:] - cell_pts[2,:]).T
- T[:,1] = (cell_pts[1,:] - cell_pts[2,:]).T
+ T = np.matrix(np.empty((2, 2)))
+ T[:, 0] = (cell_pts[0, :] - cell_pts[2, :]).T
+ T[:, 1] = (cell_pts[1, :] - cell_pts[2, :]).T
T_inv = np.linalg.inv(T)
- dl1 = T_inv[0,:] # row 0
- dl2 = T_inv[1,:] # row 1
+ dl1 = T_inv[0, :] # row 0
+ dl2 = T_inv[1, :] # row 1
for node in range(3):
l1, l2 = T_inv * (cell_pts[node, :].T - cell_pts[2, :].T)
- assert -1e-15 < l1 and 1 + 1e-15 > l1 \
- and -1e-15 < l2 and 1+ 1e-15 > l2
+ assert -1e-15 < l1 and 1 + 1e-15 > l1 and -1e-15 < l2 and 1 + 1e-15 > l2
- grad = np.empty((2,2))
+ grad = np.empty((2, 2))
for comp in range(2):
nodal_values = point_data[node_ids, comp]
# nodal_values = cell_pts[:, comp].flat
# if t > 0 and cell_pts[0,1] > 0.95 and comp == 1:
# print(nodal_values[0])
- grad[comp,:] = dl1 * nodal_values[0] \
- + dl2 * nodal_values[1] \
- - (dl1 + dl2) * nodal_values[2]
+ grad[comp, :] = (
+ dl1 * nodal_values[0]
+ + dl2 * nodal_values[1]
+ - (dl1 + dl2) * nodal_values[2]
+ )
# if t > 0 and cell_pts[0,1] > 0.95:
# print(grad)
@@ -180,15 +187,14 @@ for t, fn in zip(ts, fns):
node_id = node_ids[node]
if r == 0:
- dvdr = grad[0,0]
+ dvdr = grad[0, 0]
v_over_r = dvdr
else:
v_over_r = point_data[node_id, 0] / r
- strain_kelvin = np.array([
- strain[0,0], strain[1,1], v_over_r,
- strain[0,1] * np.sqrt(2.0)
- ])
+ strain_kelvin = np.array(
+ [strain[0, 0], strain[1, 1], v_over_r, strain[0, 1] * np.sqrt(2.0)]
+ )
strain_data[node_id, :] = strain_kelvin
def computeStrain(grid):
@@ -209,7 +215,7 @@ for t, fn in zip(ts, fns):
grid.GetPointData().AddArray(strain_kelvin_vtk)
strain = strain_kelvin.copy()
- strain[:,3] /= np.sqrt(2.0)
+ strain[:, 3] /= np.sqrt(2.0)
strain_vtk = numpy_to_vtk(strain, 1)
strain_vtk.SetName("strain_post")
grid.GetPointData().AddArray(strain_vtk)
@@ -224,15 +230,16 @@ for t, fn in zip(ts, fns):
stress_kelvin_vtk.SetName("stress_post_kelvin")
stress_symm_tensor = stress_kelvin.copy()
- stress_symm_tensor[:,3] /= np.sqrt(2.0)
+ stress_symm_tensor[:, 3] /= np.sqrt(2.0)
stress_symm_tensor_vtk = numpy_to_vtk(stress_symm_tensor, 1)
stress_symm_tensor_vtk.SetName("stress_post")
grid.GetPointData().AddArray(stress_symm_tensor_vtk)
writer.SetInputData(grid)
- writer.SetFileName(os.path.join(
- os.path.dirname(fn), "post_{:.0f}.vtu".format(t)))
+ writer.SetFileName(
+ os.path.join(os.path.dirname(fn), "post_{:.0f}.vtu".format(t))
+ )
writer.Write()
return grid
@@ -255,8 +262,8 @@ for t, fn in zip(ts, fns):
# determine top boundary
coords = vtk_to_numpy(grid.GetPoints().GetData())
- assert abs(min(coords[:,0])) < 1e-8
- idcs_top_boundary = np.where(coords[:,1] > y_top - 1e-7)[0]
+ assert abs(min(coords[:, 0])) < 1e-8
+ idcs_top_boundary = np.where(coords[:, 1] > y_top - 1e-7)[0]
# print(idcs_top_boundary)
assert len(idcs_top_boundary) != 0
xs_top_boundary = coords[idcs_top_boundary, 0]
@@ -269,14 +276,14 @@ for t, fn in zip(ts, fns):
def average_stress(rs, stress):
r_contact = max(rs)
- rs_int = np.linspace(min(rs), max(rs)-1e-8, max(len(rs), 200))
+ rs_int = np.linspace(min(rs), max(rs) - 1e-8, max(len(rs), 200))
stress_int = interp1d(rs, stress, bounds_error=False, fill_value=0.0)
avg_stress = np.empty_like(rs_int)
for i, r in enumerate(rs_int):
- rho_max = np.sqrt(r_contact**2 - r**2)
+ rho_max = np.sqrt(r_contact ** 2 - r ** 2)
rhos = np.linspace(0, rho_max, 100)
- xis = np.sqrt(rhos**2 + r**2)
+ xis = np.sqrt(rhos ** 2 + r ** 2)
try:
assert max(xis) <= r_contact + 1e-8
except:
@@ -296,7 +303,6 @@ for t, fn in zip(ts, fns):
return F
-
def stress_at_contact_area():
global add_leg
@@ -318,28 +324,41 @@ for t, fn in zip(ts, fns):
w_0 = 2.0 * (1.0 - y_top)
rs = np.linspace(0, r_contact, 200)
- if add_leg: ax.plot([], [], color="k", ls=":", label="ref")
- h, = ax.plot(rs, -p_contact(rs, r_contact), ls=":")
+ if add_leg:
+ ax.plot([], [], color="k", ls=":", label="ref")
+ (h,) = ax.plot(rs, -p_contact(rs, r_contact), ls=":")
if False:
r_contact_ana = np.sqrt(w_0 * R)
rs2 = np.linspace(0, r_contact_ana, 200)
- if add_leg: ax.plot([], [], color="k", ls="--", label="ref2")
+ if add_leg:
+ ax.plot([], [], color="k", ls="--", label="ref2")
ax.plot(rs2, -p_contact(rs, r_contact_ana), ls="--", color=h.get_color())
- stress_yy = vtk_to_numpy(grid.GetPointData().GetArray("sigma"))[:,1]
+ stress_yy = vtk_to_numpy(grid.GetPointData().GetArray("sigma"))[:, 1]
rs, avg_stress_yy = average_stress(xs, stress_yy)
- if add_leg: ax.plot([], [], color="k", ls="-", label="ogs")
- ax.plot(rs, avg_stress_yy, color=h.get_color(), ls="-",
- label=r"$w_0 = {}$".format(w_0))
+ if add_leg:
+ ax.plot([], [], color="k", ls="-", label="ogs")
+ ax.plot(
+ rs,
+ avg_stress_yy,
+ color=h.get_color(),
+ ls="-",
+ label=r"$w_0 = {}$".format(w_0),
+ )
if False:
stress = vtk_to_numpy(grid.GetPointData().GetArray("stress_post"))
- rs, avg_stress_yy = average_stress(xs, stress[:,1])
- if add_leg: ax.plot([], [], color="k", label="post")
- ax.plot(rs, avg_stress_yy, color=h.get_color(),
- label=r"$w_0 = {}$".format(2*(1.0 - y_top)))
+ rs, avg_stress_yy = average_stress(xs, stress[:, 1])
+ if add_leg:
+ ax.plot([], [], color="k", label="post")
+ ax.plot(
+ rs,
+ avg_stress_yy,
+ color=h.get_color(),
+ label=r"$w_0 = {}$".format(2 * (1.0 - y_top)),
+ )
ax.scatter([r_contact], [0], color=h.get_color())
ax.legend(loc="center left", bbox_to_anchor=(1.0, 0.5))
@@ -378,12 +397,12 @@ fig, ax = plt.subplots()
ax.scatter(rs_contact[1:], Fs, label="ogs")
rs_ref = np.linspace(0, max(rs_contact), 200)
-Fs_ref = 8. * rs_ref**3 * kappa / 3.0
+Fs_ref = 8.0 * rs_ref ** 3 * kappa / 3.0
ax.plot(rs_ref, Fs_ref, label="ref")
l = ax.legend()
-l.get_frame().set_facecolor('white')
+l.get_frame().set_facecolor("white")
ax.set_xlabel("radius of contact area $a$ / m")
ax.set_ylabel("applied force $F$ / N")
diff --git a/Tests/Data/Mechanics/Linear/PythonPiston/chamber.py b/Tests/Data/Mechanics/Linear/PythonPiston/chamber.py
index 2fd26893826..d5639d42598 100644
--- a/Tests/Data/Mechanics/Linear/PythonPiston/chamber.py
+++ b/Tests/Data/Mechanics/Linear/PythonPiston/chamber.py
@@ -5,15 +5,21 @@ from math import pi
R_piston = 0.05
h_chamber_0 = 0.1
-V_chamber_0 = pi * R_piston**2 * h_chamber_0
+V_chamber_0 = pi * R_piston ** 2 * h_chamber_0
p_chamber_0 = 1e5
def p_chamber(displacement_y):
- V_chamber = pi * R_piston**2 * (h_chamber_0 + displacement_y)
+ V_chamber = pi * R_piston ** 2 * (h_chamber_0 + displacement_y)
return p_chamber_0 * V_chamber_0 / V_chamber
+
def dp_chamber_du_y(displacement_y):
- return -p_chamber_0 * V_chamber_0 / pi / R_piston**2 / \
- (h_chamber_0 + displacement_y)**2
+ return (
+ -p_chamber_0
+ * V_chamber_0
+ / pi
+ / R_piston ** 2
+ / (h_chamber_0 + displacement_y) ** 2
+ )
diff --git a/Tests/Data/Mechanics/Linear/PythonPiston/piston_bc.py b/Tests/Data/Mechanics/Linear/PythonPiston/piston_bc.py
index e06fbf50238..b7437e06bc5 100644
--- a/Tests/Data/Mechanics/Linear/PythonPiston/piston_bc.py
+++ b/Tests/Data/Mechanics/Linear/PythonPiston/piston_bc.py
@@ -23,7 +23,6 @@ class BC_Y(OpenGeoSys.BoundaryCondition):
return (False, 0.0)
-
def getFlux(self, t, coords, primary_vars):
x, y, z = coords
@@ -32,10 +31,10 @@ class BC_Y(OpenGeoSys.BoundaryCondition):
ux, uy = primary_vars
p = p_chamber(uy)
- Jac = [ 0.0, dp_chamber_du_y(uy) ]
+ Jac = [0.0, dp_chamber_du_y(uy)]
return (True, p, Jac)
- return (False, 0.0, [ 0.0, 0.0 ])
+ return (False, 0.0, [0.0, 0.0])
# instantiate the BC objects used by OpenGeoSys
diff --git a/Tests/Data/Mechanics/Linear/PythonPiston/post.py b/Tests/Data/Mechanics/Linear/PythonPiston/post.py
index 09695a52a3b..73a63e93685 100755
--- a/Tests/Data/Mechanics/Linear/PythonPiston/post.py
+++ b/Tests/Data/Mechanics/Linear/PythonPiston/post.py
@@ -15,6 +15,7 @@ pvd_file = "out/piston_pcs_0.pvd"
### helpers ##############################################
import os
+
try:
import xml.etree.cElementTree as ET
except:
@@ -26,6 +27,7 @@ def relpathfrom(origin, relpath):
return relpath
return os.path.join(origin, relpath)
+
def read_pvd_file(fn):
try:
path = fn.name
@@ -34,22 +36,27 @@ def read_pvd_file(fn):
pathroot = os.path.dirname(path)
pvdtree = ET.parse(fn)
node = pvdtree.getroot()
- if node.tag != "VTKFile": return None, None
+ if node.tag != "VTKFile":
+ return None, None
children = list(node)
- if len(children) != 1: return None, None
+ if len(children) != 1:
+ return None, None
node = children[0]
- if node.tag != "Collection": return None, None
+ if node.tag != "Collection":
+ return None, None
ts = []
fs = []
for child in node:
- if child.tag != "DataSet": return None, None
+ if child.tag != "DataSet":
+ return None, None
ts.append(float(child.get("timestep")))
fs.append(relpathfrom(pathroot, child.get("file")))
return ts, fs
+
### helpers end ##########################################
@@ -59,7 +66,7 @@ reader = vtk.vtkXMLUnstructuredGridReader()
loc_points = vtk.vtkPoints()
-loc_points.InsertNextPoint([0.0, 0.0, 0.0])
+loc_points.InsertNextPoint([0.0, 0.0, 0.0])
loc = vtk.vtkPolyData()
loc.SetPoints(loc_points)
@@ -76,8 +83,8 @@ for i, (t, fn) in enumerate(zip(ts, fns)):
probe.Update()
grid = probe.GetOutput()
- uy = vtk_to_numpy(grid.GetPointData().GetArray("displacement"))[0,1]
- p = vtk_to_numpy(grid.GetPointData().GetArray("sigma"))[0,1]
+ uy = vtk_to_numpy(grid.GetPointData().GetArray("displacement"))[0, 1]
+ p = vtk_to_numpy(grid.GetPointData().GetArray("sigma"))[0, 1]
uys[i] = uy
ps[i] = -p
diff --git a/Tests/Data/Parabolic/T/1D_neumann/plotLine.py b/Tests/Data/Parabolic/T/1D_neumann/plotLine.py
index ff8f3c67fa7..4eafb6cdab4 100755
--- a/Tests/Data/Parabolic/T/1D_neumann/plotLine.py
+++ b/Tests/Data/Parabolic/T/1D_neumann/plotLine.py
@@ -35,84 +35,90 @@ def convertPointDataToPandasDF(line_data):
data[k] = vtk_to_numpy(v)
coords = line_data.GetPoints()[:, 0] # x-coordinates
- df = pd.concat([pd.DataFrame(v) for k, v in data.items()],
- axis=1,
- keys=list(data.keys()))
+ df = pd.concat(
+ [pd.DataFrame(v) for k, v in data.items()], axis=1, keys=list(data.keys())
+ )
df.index = coords
return df
-def setupMatplotlib(height=8., width=6.):
- plt.rcParams['xtick.direction'] = 'out'
- plt.rcParams['ytick.direction'] = 'out'
- plt.rcParams['lines.linewidth'] = 2.0
- plt.rcParams['lines.color'] = 'black'
- plt.rcParams['legend.frameon'] = True
- plt.rcParams['font.family'] = 'serif'
- plt.rcParams['legend.fontsize'] = 8
- plt.rcParams['font.size'] = 12
+def setupMatplotlib(height=8.0, width=6.0):
+ plt.rcParams["xtick.direction"] = "out"
+ plt.rcParams["ytick.direction"] = "out"
+ plt.rcParams["lines.linewidth"] = 2.0
+ plt.rcParams["lines.color"] = "black"
+ plt.rcParams["legend.frameon"] = True
+ plt.rcParams["font.family"] = "serif"
+ plt.rcParams["legend.fontsize"] = 8
+ plt.rcParams["font.size"] = 12
# For ipython notebook display set default values.
- #plt.rcParams['lines.markersize'] = 12
- plt.rcParams['figure.figsize'] = (height, width)
- plt.rcParams['grid.linewidth'] = 1
+ # plt.rcParams['lines.markersize'] = 12
+ plt.rcParams["figure.figsize"] = (height, width)
+ plt.rcParams["grid.linewidth"] = 1
# General settings used by display and print contexts.
- plt.rcParams['axes.axisbelow'] = True
- grid_line_color = '0.95'
- plt.rcParams['grid.color'] = grid_line_color
- plt.rcParams['grid.linestyle'] = '-'
- plt.rc('text', usetex=True)
- plt.rc('font', family='serif')
+ plt.rcParams["axes.axisbelow"] = True
+ grid_line_color = "0.95"
+ plt.rcParams["grid.color"] = grid_line_color
+ plt.rcParams["grid.linestyle"] = "-"
+ plt.rc("text", usetex=True)
+ plt.rc("font", family="serif")
def commonFormat(ax_el, centerx=None, centery=None):
ax_el.grid(True)
- ax_el.spines['top'].set_visible(0)
- ax_el.spines['right'].set_visible(0)
- ax_el.spines['bottom'].set_linewidth(0.5)
- ax_el.spines['left'].set_linewidth(0.5)
- ax_el.xaxis.set_ticks_position('bottom')
- ax_el.yaxis.set_ticks_position('left')
- if ((centerx is not None) and (centery is not None)):
- ax_el.spines['left'].set_position(('data', centerx))
- ax_el.spines['bottom'].set_position(('data', centery))
- ax_el.spines['right'].set_position(('data', centerx - 1))
- ax_el.spines['top'].set_position(('data', centery - 1))
+ ax_el.spines["top"].set_visible(0)
+ ax_el.spines["right"].set_visible(0)
+ ax_el.spines["bottom"].set_linewidth(0.5)
+ ax_el.spines["left"].set_linewidth(0.5)
+ ax_el.xaxis.set_ticks_position("bottom")
+ ax_el.yaxis.set_ticks_position("left")
+ if (centerx is not None) and (centery is not None):
+ ax_el.spines["left"].set_position(("data", centerx))
+ ax_el.spines["bottom"].set_position(("data", centery))
+ ax_el.spines["right"].set_position(("data", centerx - 1))
+ ax_el.spines["top"].set_position(("data", centery - 1))
def plotWithFormatting(ax, data, formatting, with_labels=True):
- values = data[formatting['column']]
- if 'component' in formatting:
- values = values[formatting['component']]
+ values = data[formatting["column"]]
+ if "component" in formatting:
+ values = values[formatting["component"]]
- #print("plot", formatting['column'], values)
- return ax.plot(data.index,
- values,
- ls=formatting['ls'],
- color=formatting['color'],
- lw=1.,
- label=formatting['label'] if with_labels else '')
+ # print("plot", formatting['column'], values)
+ return ax.plot(
+ data.index,
+ values,
+ ls=formatting["ls"],
+ color=formatting["color"],
+ lw=1.0,
+ label=formatting["label"] if with_labels else "",
+ )
def plotTemperature(ax, data, with_labels=True):
- formattings = [{
- 'column': ('newton', 'temperature'),
- 'component': 0,
- 'color': 'green',
- 'ls': '-',
- 'label': 'no mass lumping'
- }, {
- 'column': ('newton_masslumping', 'temperature'),
- 'component': 0,
- 'color': 'blue',
- 'ls': '-',
- 'label': 'with mass lumping'
- }, {
- 'column': ('newton', 'analytical'),
- 'color': '#e41a1c',
- 'ls': '-',
- 'label': 'analytical'
- }]
+ formattings = [
+ {
+ "column": ("newton", "temperature"),
+ "component": 0,
+ "color": "green",
+ "ls": "-",
+ "label": "no mass lumping",
+ },
+ {
+ "column": ("newton_masslumping", "temperature"),
+ "component": 0,
+ "color": "blue",
+ "ls": "-",
+ "label": "with mass lumping",
+ },
+ {
+ "column": ("newton", "analytical"),
+ "color": "#e41a1c",
+ "ls": "-",
+ "label": "analytical",
+ },
+ ]
lines = []
for formatting in formattings:
lines += plotWithFormatting(ax, data, formatting, with_labels)
@@ -120,17 +126,20 @@ def plotTemperature(ax, data, with_labels=True):
def plotErrors(ax, data, with_labels=True):
- formattings = [{
- 'column': ('newton', 'error'),
- 'color': 'green',
- 'ls': '-',
- 'label': '$e_\mathrm{newton}$'
- }, {
- 'column': ('newton_masslumping', 'error'),
- 'color': 'blue',
- 'ls': '-',
- 'label': '$e_\mathrm{newton}^\mathrm{ML}$'
- }]
+ formattings = [
+ {
+ "column": ("newton", "error"),
+ "color": "green",
+ "ls": "-",
+ "label": "$e_\mathrm{newton}$",
+ },
+ {
+ "column": ("newton_masslumping", "error"),
+ "color": "blue",
+ "ls": "-",
+ "label": "$e_\mathrm{newton}^\mathrm{ML}$",
+ },
+ ]
lines = []
for formatting in formattings:
lines += plotWithFormatting(ax, data, formatting, with_labels)
@@ -138,46 +147,46 @@ def plotErrors(ax, data, with_labels=True):
def plotCases(data, output):
- d = data[data[('picard_masslumping', 'error')].abs() > 1e-6]
+ d = data[data[("picard_masslumping", "error")].abs() > 1e-6]
setupMatplotlib(6, 7)
- plt.rcParams['legend.fontsize'] = 10
- plt.rcParams['font.size'] = 14
+ plt.rcParams["legend.fontsize"] = 10
+ plt.rcParams["font.size"] = 14
fig, axes = plt.subplots(nrows=2, ncols=1)
lines = plotTemperature(axes[0], d, with_labels=True)
- axes[0].set_ylabel('Temperature / K', fontsize=18)
- axes[0].legend(loc='best',
- fontsize=12,
- ncol=1,
- bbox_to_anchor=(0.5, 0.5, 0.5, 0.5))
+ axes[0].set_ylabel("Temperature / K", fontsize=18)
+ axes[0].legend(loc="best", fontsize=12, ncol=1, bbox_to_anchor=(0.5, 0.5, 0.5, 0.5))
lines += plotErrors(axes[1], d, with_labels=False)
- axes[1].set_ylabel('Error / K', fontsize=18)
- axes[1].set_xlabel('$x$ / m', fontsize=16)
+ axes[1].set_ylabel("Error / K", fontsize=18)
+ axes[1].set_xlabel("$x$ / m", fontsize=16)
for ax in axes:
commonFormat(ax)
- #ax.set_xlim(left=0.)
- #ax.set_xlim(right=5.)
+ # ax.set_xlim(left=0.)
+ # ax.set_xlim(right=5.)
fig.tight_layout()
- fig.savefig(output + '.png', dpi=150)
- plt.close('all')
+ fig.savefig(output + ".png", dpi=150)
+ plt.close("all")
return None
def plotPicardTemperature(ax, data, with_labels=True):
- formattings = [{
- 'column': ('picard', 'temperature'),
- 'component': 0,
- 'color': 'blue',
- 'ls': '-',
- 'label': 'picard, no mass lumping'
- }, {
- 'column': ('picard_masslumping', 'temperature'),
- 'component': 0,
- 'color': 'green',
- 'ls': '-',
- 'label': 'picard, with mass lumping'
- }]
+ formattings = [
+ {
+ "column": ("picard", "temperature"),
+ "component": 0,
+ "color": "blue",
+ "ls": "-",
+ "label": "picard, no mass lumping",
+ },
+ {
+ "column": ("picard_masslumping", "temperature"),
+ "component": 0,
+ "color": "green",
+ "ls": "-",
+ "label": "picard, with mass lumping",
+ },
+ ]
lines = []
for formatting in formattings:
lines += plotWithFormatting(ax, data, formatting, with_labels)
@@ -185,19 +194,22 @@ def plotPicardTemperature(ax, data, with_labels=True):
def plotNewtonTemperature(ax, data, with_labels=True):
- formattings = [{
- 'column': ('newton', 'temperature'),
- 'component': 0,
- 'color': 'orange',
- 'ls': '-.',
- 'label': 'newton, no mass lumping'
- }, {
- 'column': ('newton_masslumping', 'temperature'),
- 'component': 0,
- 'color': 'red',
- 'ls': '-.',
- 'label': 'newton, with mass lumping'
- }]
+ formattings = [
+ {
+ "column": ("newton", "temperature"),
+ "component": 0,
+ "color": "orange",
+ "ls": "-.",
+ "label": "newton, no mass lumping",
+ },
+ {
+ "column": ("newton_masslumping", "temperature"),
+ "component": 0,
+ "color": "red",
+ "ls": "-.",
+ "label": "newton, with mass lumping",
+ },
+ ]
lines = []
for formatting in formattings:
lines += plotWithFormatting(ax, data, formatting, with_labels)
@@ -205,17 +217,10 @@ def plotNewtonTemperature(ax, data, with_labels=True):
def plotNewtonVsPicardErrors(ax, data, with_labels=True):
- formattings = [{
- 'column': 'no_ml',
- 'color': 'green',
- 'ls': '-',
- 'label': 'no mass lumping'
- }, {
- 'column': 'ml',
- 'color': 'blue',
- 'ls': '-',
- 'label': 'with mass lumping'
- }]
+ formattings = [
+ {"column": "no_ml", "color": "green", "ls": "-", "label": "no mass lumping"},
+ {"column": "ml", "color": "blue", "ls": "-", "label": "with mass lumping"},
+ ]
lines = []
for formatting in formattings:
lines += plotWithFormatting(ax, data, formatting, with_labels)
@@ -223,62 +228,66 @@ def plotNewtonVsPicardErrors(ax, data, with_labels=True):
def plotNewtonVsPicard(data, output):
- data['no_ml'] = data[('newton', 'temperature',
- 0)] - data[('picard', 'temperature', 0)]
- data['ml'] = data[('newton_masslumping', 'temperature',
- 0)] - data[('picard_masslumping', 'temperature', 0)]
- #data['select'] = (data['ml'].abs() > 1e-16) | (data['no_ml'].abs() > 1e-16)
- #print(d[d].index[-1])
+ data["no_ml"] = (
+ data[("newton", "temperature", 0)] - data[("picard", "temperature", 0)]
+ )
+ data["ml"] = (
+ data[("newton_masslumping", "temperature", 0)]
+ - data[("picard_masslumping", "temperature", 0)]
+ )
+ # data['select'] = (data['ml'].abs() > 1e-16) | (data['no_ml'].abs() > 1e-16)
+ # print(d[d].index[-1])
setupMatplotlib(6, 7)
- plt.rcParams['legend.fontsize'] = 10
- plt.rcParams['font.size'] = 14
+ plt.rcParams["legend.fontsize"] = 10
+ plt.rcParams["font.size"] = 14
fig, axes = plt.subplots(nrows=2, ncols=1)
lines = plotPicardTemperature(axes[0], data, with_labels=True)
lines = plotNewtonTemperature(axes[0], data, with_labels=True)
- axes[0].set_ylabel('Temperature / K', fontsize=18)
- axes[0].legend(loc='best',
- fontsize=12,
- ncol=1,
- bbox_to_anchor=(0.5, 0.5, 0.5, 0.5))
+ axes[0].set_ylabel("Temperature / K", fontsize=18)
+ axes[0].legend(loc="best", fontsize=12, ncol=1, bbox_to_anchor=(0.5, 0.5, 0.5, 0.5))
lines += plotNewtonVsPicardErrors(axes[1], data, with_labels=True)
- axes[1].legend(loc='best',
- fontsize=12,
- ncol=1,
- bbox_to_anchor=(0.5, 0.5, 0.5, 0.5))
- axes[1].set_ylabel('$T_\mathrm{diff}$ (newton - picard) / K', fontsize=18)
- axes[1].set_xlabel('$x$ / m', fontsize=16)
+ axes[1].legend(loc="best", fontsize=12, ncol=1, bbox_to_anchor=(0.5, 0.5, 0.5, 0.5))
+ axes[1].set_ylabel("$T_\mathrm{diff}$ (newton - picard) / K", fontsize=18)
+ axes[1].set_xlabel("$x$ / m", fontsize=16)
for ax in axes:
commonFormat(ax)
- #ax.set_xlim(left=0.)
- #ax.set_xlim(right=5.)
+ # ax.set_xlim(left=0.)
+ # ax.set_xlim(right=5.)
fig.tight_layout()
- fig.savefig(output + '.png', dpi=300)
- plt.close('all')
+ fig.savefig(output + ".png", dpi=300)
+ plt.close("all")
return None
def singleTimeStep(ts):
dt = 78125
t = ts * dt
- cases = ['picard', 'newton', 'picard_masslumping', 'newton_masslumping']
+ cases = ["picard", "newton", "picard_masslumping", "newton_masslumping"]
data = []
for c in cases:
- line = probeFileAlongLine(c + '_ts_' + str(ts) + '_t_' + str(t) +
- '.000000.vtu')
+ line = probeFileAlongLine(c + "_ts_" + str(ts) + "_t_" + str(t) + ".000000.vtu")
df = convertPointDataToPandasDF(line)
- df['analytical'] = temperature_analytical.temperature(df.index, t)
- df['error'] = df['temperature'].sub(df['analytical'], axis=0)
- print('Error for ts', ts, 'case', c, 'is',
- np.sqrt((df['error']**2).sum() * 60 / 10000))
- data.append(df[['temperature', 'error', 'analytical']])
-
- df = pd.concat(data, keys=cases, names=['case'], axis=1)
- plotCases(df[df[('picard_masslumping', 'error')].abs() > 1e-6],
- 'temperature_error_ts_' + str(ts) + '_t_' + str(t))
- plotNewtonVsPicard(df, 'picard_vs_newton_ts_' + str(ts) + '_t_' + str(t))
+ df["analytical"] = temperature_analytical.temperature(df.index, t)
+ df["error"] = df["temperature"].sub(df["analytical"], axis=0)
+ print(
+ "Error for ts",
+ ts,
+ "case",
+ c,
+ "is",
+ np.sqrt((df["error"] ** 2).sum() * 60 / 10000),
+ )
+ data.append(df[["temperature", "error", "analytical"]])
+
+ df = pd.concat(data, keys=cases, names=["case"], axis=1)
+ plotCases(
+ df[df[("picard_masslumping", "error")].abs() > 1e-6],
+ "temperature_error_ts_" + str(ts) + "_t_" + str(t),
+ )
+ plotNewtonVsPicard(df, "picard_vs_newton_ts_" + str(ts) + "_t_" + str(t))
def main():
diff --git a/Tests/Data/Parabolic/T/1D_neumann/temperature_analytical.py b/Tests/Data/Parabolic/T/1D_neumann/temperature_analytical.py
index 1e09f0c3320..1aa3a3eabf6 100755
--- a/Tests/Data/Parabolic/T/1D_neumann/temperature_analytical.py
+++ b/Tests/Data/Parabolic/T/1D_neumann/temperature_analytical.py
@@ -6,7 +6,7 @@ import numpy as np
from scipy.special import erfc
r = vtkXMLUnstructuredGridReader()
-r.SetFileName('mesh.vtu')
+r.SetFileName("mesh.vtu")
r.Update()
m = r.GetOutput()
ps = m.GetPoints()
@@ -23,13 +23,15 @@ q = 2
def temperature(x, t):
- return T_inf + 2 * q / lambda_coeff * (np.sqrt(alpha * t / np.pi) * np.exp(
- -x**2 / (4 * alpha * t)) - x / 2 * erfc(x / (2 * np.sqrt(alpha * t))))
+ return T_inf + 2 * q / lambda_coeff * (
+ np.sqrt(alpha * t / np.pi) * np.exp(-(x ** 2) / (4 * alpha * t))
+ - x / 2 * erfc(x / (2 * np.sqrt(alpha * t)))
+ )
def addSolution(t):
T = vtkDoubleArray()
- T.SetName('temperature_' + str(t) + 's')
+ T.SetName("temperature_" + str(t) + "s")
T.SetNumberOfComponents(1)
T.SetNumberOfTuples(ps.GetNumberOfPoints())
@@ -44,6 +46,6 @@ for t in 78125 * np.array([1, 3, 65, 405, 500]):
addSolution(t)
w = vtkXMLUnstructuredGridWriter()
-w.SetFileName('x.vtu')
+w.SetFileName("x.vtu")
w.SetInputData(m)
w.Update()
diff --git a/Tests/Data/Parabolic/T/3D_3BHEs_array/bcs_tespy.py b/Tests/Data/Parabolic/T/3D_3BHEs_array/bcs_tespy.py
index 9451cf8fc35..b9d1f6f1aa5 100644
--- a/Tests/Data/Parabolic/T/3D_3BHEs_array/bcs_tespy.py
+++ b/Tests/Data/Parabolic/T/3D_3BHEs_array/bcs_tespy.py
@@ -6,6 +6,7 @@
###
import sys
+
print(sys.version)
import os
import numpy as np
@@ -19,14 +20,14 @@ from tespy.networks import load_network
rho_f = 992.92 # kg/m3
# switch for special boundary conditions
# switch of the function for manually specified dynamic flowrate
-switch_dyn_frate = 'off' # 'on','off'
+switch_dyn_frate = "off" # 'on','off'
# switch of the function for manually specified dynamic thermal demand
-switch_dyn_demand = 'on' # 'on','off'
+switch_dyn_demand = "on" # 'on','off'
# network status setting
def network_status(t):
- nw_status = 'on'
+ nw_status = "on"
# month for closed network
timerange_nw_off_month = [] # No month for closed network
# t-1 to avoid the calculation problem at special time point,
@@ -36,7 +37,7 @@ def network_status(t):
if t_trans_month > 12:
t_trans_month = t_trans - 12 * (int(t_trans / 12))
if t_trans_month in timerange_nw_off_month:
- nw_status = 'off'
+ nw_status = "off"
return nw_status
@@ -49,9 +50,18 @@ def consumer_demand(t): # dynamic thermal demand from consumer
# thermal demand in each month (assumed specific heat extraction rate*
# length of BHE* number of BHE)
month_demand = [
- -25 * 50 * 3, -25 * 50 * 3, -25 * 50 * 3, -25 * 50 * 3, -25 * 50 * 3,
- -25 * 50 * 3, -25 * 50 * 3, -25 * 50 * 3, -25 * 50 * 3, -25 * 50 * 3,
- -25 * 50 * 3, -25 * 50 * 3
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
]
return month_demand[t_trans - 1]
@@ -73,42 +83,42 @@ def dyn_frate(t):
# create network dataframe
def create_dataframe():
# return dataframe
- df_nw = read_csv('./pre/bhe_network.csv',
- delimiter=';',
- index_col=[0],
- dtype={'data_index': str})
- return (df_nw)
+ df_nw = read_csv(
+ "./pre/bhe_network.csv", delimiter=";", index_col=[0], dtype={"data_index": str}
+ )
+ return df_nw
# TESPy calculation process
def get_tespy_results(t):
# bhe network boundary conditions re parametrization
# if network exist dynamic flowrate
- if switch_dyn_frate == 'on':
+ if switch_dyn_frate == "on":
cur_frate = dyn_frate(t)
- localVars['inlet_name'].set_attr(m=cur_frate)
+ localVars["inlet_name"].set_attr(m=cur_frate)
# if network exist dynamic thermal demand
- if switch_dyn_demand == 'on':
+ if switch_dyn_demand == "on":
# consumer thermal load:
cur_month_demand = consumer_demand(t)
nw.busses[bus_name].set_attr(P=cur_month_demand)
# T_out re parametrization:
for i in range(n_BHE):
- localVars['outlet_BHE' + str(i + 1)].set_attr(T=df.loc[data_index[i],
- 'Tout_val'])
+ localVars["outlet_BHE" + str(i + 1)].set_attr(
+ T=df.loc[data_index[i], "Tout_val"]
+ )
# solving network
- nw.solve(mode='design')
+ nw.solve(mode="design")
# get Tin_val and flow rate
for i in range(n_BHE):
# get Tin_val
- df.loc[df.index[i],
- 'Tin_val'] = localVars['inlet_BHE' +
- str(i + 1)].get_attr('T').val
+ df.loc[df.index[i], "Tin_val"] = (
+ localVars["inlet_BHE" + str(i + 1)].get_attr("T").val
+ )
# get flowrate
- df.loc[df.index[i],
- 'flowrate'] = localVars['inlet_BHE' +
- str(i + 1)].get_attr('m').val / rho_f
- return df['Tin_val'].tolist(), df['flowrate'].tolist()
+ df.loc[df.index[i], "flowrate"] = (
+ localVars["inlet_BHE" + str(i + 1)].get_attr("m").val / rho_f
+ )
+ return df["Tin_val"].tolist(), df["flowrate"].tolist()
# OGS setting
@@ -116,35 +126,36 @@ def get_tespy_results(t):
class BC(OpenGeoSys.BHENetwork):
def initializeDataContainer(self):
# initialize network and get data from the network
- nw.solve(mode='design')
+ nw.solve(mode="design")
get_tespy_results(0)
# convert dataframe to column list
t = 0 # 'initial time'
- data_col_1 = df['Tin_val'].tolist() # 'Tin_val'
- data_col_2 = df['Tout_val'].tolist() # 'Tout_val'
- data_col_3 = df['Tout_node_id'].astype(int).tolist() # 'Tout_node_id'
- data_col_4 = df['flowrate'].tolist() # 'BHE flow rate'
+ data_col_1 = df["Tin_val"].tolist() # 'Tin_val'
+ data_col_2 = df["Tout_val"].tolist() # 'Tout_val'
+ data_col_3 = df["Tout_node_id"].astype(int).tolist() # 'Tout_node_id'
+ data_col_4 = df["flowrate"].tolist() # 'BHE flow rate'
return (t, data_col_1, data_col_2, data_col_3, data_col_4)
def tespySolver(self, t, Tin_val, Tout_val):
# network status:
nw_status = network_status(t)
# if network closed:
- if nw_status == 'off':
- df.loc[:,'flowrate'] = 0
- cur_flowrate = df['flowrate'].tolist()
+ if nw_status == "off":
+ df.loc[:, "flowrate"] = 0
+ cur_flowrate = df["flowrate"].tolist()
return (True, True, Tout_val, cur_flowrate)
else:
# read Tout_val to dataframe
for i in range(n_BHE):
- df.loc[df.index[i], 'Tout_val'] = Tout_val[i]
+ df.loc[df.index[i], "Tout_val"] = Tout_val[i]
# TESPy solver
cur_Tin_val, cur_flowrate = get_tespy_results(t)
# check norm if network achieves the converge
if_success = False
pre_Tin_val = Tin_val
norm = np.linalg.norm(
- abs(np.asarray(pre_Tin_val) - np.asarray(cur_Tin_val)))
+ abs(np.asarray(pre_Tin_val) - np.asarray(cur_Tin_val))
+ )
if norm < 10e-6:
if_success = True
# return to OGS
@@ -156,7 +167,7 @@ class BC(OpenGeoSys.BHENetwork):
# load path of network model:
# loading the TESPy model
project_dir = os.chdir(ogs_prj_directory)
-nw = load_network('./pre/tespy_nw')
+nw = load_network("./pre/tespy_nw")
# set if print the network iteration info
nw.set_attr(iterinfo=False)
@@ -172,28 +183,28 @@ for i in range(n_BHE):
for c in nw.conns.index:
# bhe inlet and outlet conns
if c.target.label == data_index[i]: # inlet conns of bhe
- localVars['inlet_BHE' + str(i + 1)] = c
+ localVars["inlet_BHE" + str(i + 1)] = c
if c.source.label == data_index[i]: # outlet conns of bhe
- localVars['outlet_BHE' + str(i + 1)] = c
+ localVars["outlet_BHE" + str(i + 1)] = c
# time depended flowrate
-if switch_dyn_frate == 'on':
+if switch_dyn_frate == "on":
# import the name of inlet connection from the network csv file
- inlet_name = read_csv('./pre/tespy_nw/connections.csv',
- delimiter=';',
- index_col=[0]).iloc[0,0]
+ inlet_name = read_csv(
+ "./pre/tespy_nw/connections.csv", delimiter=";", index_col=[0]
+ ).iloc[0, 0]
for c in nw.conns.index:
# bhe inflow conns
if c.source.label == inlet_name: # inlet conns of bhe
- localVars['inlet_name'] = c
+ localVars["inlet_name"] = c
# time depended consumer thermal demand
-if switch_dyn_demand == 'on':
+if switch_dyn_demand == "on":
# import the name of bus from the network csv file
- bus_name = read_csv('./pre/tespy_nw/components/bus.csv',
- delimiter=';',
- index_col=[0]).index[0]
+ bus_name = read_csv(
+ "./pre/tespy_nw/components/bus.csv", delimiter=";", index_col=[0]
+ ).index[0]
# instantiate BC objects referenced in OpenGeoSys
bc_bhe = BC()
diff --git a/Tests/Data/Parabolic/T/3D_3BHEs_array/bcs_tespy_closedloop.py b/Tests/Data/Parabolic/T/3D_3BHEs_array/bcs_tespy_closedloop.py
index dd975c35043..8765c0b7639 100644
--- a/Tests/Data/Parabolic/T/3D_3BHEs_array/bcs_tespy_closedloop.py
+++ b/Tests/Data/Parabolic/T/3D_3BHEs_array/bcs_tespy_closedloop.py
@@ -6,6 +6,7 @@
###
import sys
+
print(sys.version)
import os
import numpy as np
@@ -19,16 +20,16 @@ from tespy.networks import load_network
rho_f = 992.92 # kg/m3
# switch for special boundary conditions
# switch of the function for manually specified dynamic flowrate
-switch_dyn_frate = 'off' # 'on','off'
+switch_dyn_frate = "off" # 'on','off'
# switch of the function for manually specified dynamic thermal demand
-switch_dyn_demand = 'on' # 'on','off'
-if switch_dyn_demand == 'on':
- #give the consumer name defined by user in the network model
- consumer_name = 'consumer'
+switch_dyn_demand = "on" # 'on','off'
+if switch_dyn_demand == "on":
+ # give the consumer name defined by user in the network model
+ consumer_name = "consumer"
# network status setting
def network_status(t):
- nw_status = 'on'
+ nw_status = "on"
# month for closed network
timerange_nw_off_month = [] # No month for closed network
# t-1 to avoid the calculation problem at special time point,
@@ -38,7 +39,7 @@ def network_status(t):
if t_trans_month > 12:
t_trans_month = t_trans - 12 * (int(t_trans / 12))
if t_trans_month in timerange_nw_off_month:
- nw_status = 'off'
+ nw_status = "off"
return nw_status
@@ -51,9 +52,18 @@ def consumer_demand(t): # dynamic thermal demand from consumer
# thermal demand in each month (assumed specific heat extraction rate*
# length of BHE* number of BHE)
month_demand = [
- -25 * 50 * 3, -25 * 50 * 3, -25 * 50 * 3, -25 * 50 * 3, -25 * 50 * 3,
- -25 * 50 * 3, -25 * 50 * 3, -25 * 50 * 3, -25 * 50 * 3, -25 * 50 * 3,
- -25 * 50 * 3, -25 * 50 * 3
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
+ -25 * 50 * 3,
]
return month_demand[t_trans - 1]
@@ -75,42 +85,42 @@ def dyn_frate(t):
# create network dataframe
def create_dataframe():
# return dataframe
- df_nw = read_csv('./pre/bhe_network.csv',
- delimiter=';',
- index_col=[0],
- dtype={'data_index': str})
- return (df_nw)
+ df_nw = read_csv(
+ "./pre/bhe_network.csv", delimiter=";", index_col=[0], dtype={"data_index": str}
+ )
+ return df_nw
# TESPy calculation process
def get_tespy_results(t):
# bhe network boundary conditions re parametrization
# if network exist dynamic flowrate
- if switch_dyn_frate == 'on':
+ if switch_dyn_frate == "on":
cur_frate = dyn_frate(t)
- localVars['inlet_name'].set_attr(m=cur_frate)
+ localVars["inlet_name"].set_attr(m=cur_frate)
# if network exist dynamic thermal demand
- if switch_dyn_demand == 'on':
+ if switch_dyn_demand == "on":
# consumer thermal load:
cur_month_demand = consumer_demand(t)
nw.components[consumer_name].set_attr(Q=cur_month_demand)
# T_out re parametrization:
for i in range(n_BHE):
- localVars['outlet_BHE' + str(i + 1)].set_attr(T=df.loc[data_index[i],
- 'Tout_val'])
+ localVars["outlet_BHE" + str(i + 1)].set_attr(
+ T=df.loc[data_index[i], "Tout_val"]
+ )
# solving network
- nw.solve(mode='design')
+ nw.solve(mode="design")
# get Tin_val and flow rate
for i in range(n_BHE):
# get Tin_val
- df.loc[df.index[i],
- 'Tin_val'] = localVars['inlet_BHE' +
- str(i + 1)].get_attr('T').val
+ df.loc[df.index[i], "Tin_val"] = (
+ localVars["inlet_BHE" + str(i + 1)].get_attr("T").val
+ )
# get flowrate
- df.loc[df.index[i],
- 'flowrate'] = localVars['inlet_BHE' +
- str(i + 1)].get_attr('m').val / rho_f
- return df['Tin_val'].tolist(), df['flowrate'].tolist()
+ df.loc[df.index[i], "flowrate"] = (
+ localVars["inlet_BHE" + str(i + 1)].get_attr("m").val / rho_f
+ )
+ return df["Tin_val"].tolist(), df["flowrate"].tolist()
# OGS setting
@@ -118,37 +128,38 @@ def get_tespy_results(t):
class BC(OpenGeoSys.BHENetwork):
def initializeDataContainer(self):
# initialize network and get data from the network
- nw.solve(mode='design')
+ nw.solve(mode="design")
get_tespy_results(0)
# convert dataframe to column list
t = 0 # 'initial time'
- data_col_1 = df['Tin_val'].tolist() # 'Tin_val'
- data_col_2 = df['Tout_val'].tolist() # 'Tout_val'
- data_col_3 = df['Tout_node_id'].astype(int).tolist() # 'Tout_node_id'
- data_col_4 = df['flowrate'].tolist() # 'BHE flow rate'
+ data_col_1 = df["Tin_val"].tolist() # 'Tin_val'
+ data_col_2 = df["Tout_val"].tolist() # 'Tout_val'
+ data_col_3 = df["Tout_node_id"].astype(int).tolist() # 'Tout_node_id'
+ data_col_4 = df["flowrate"].tolist() # 'BHE flow rate'
return (t, data_col_1, data_col_2, data_col_3, data_col_4)
def tespySolver(self, t, Tin_val, Tout_val):
# network status:
nw_status = network_status(t)
# if network closed:
- if nw_status == 'off':
- df.loc[:,'flowrate'] = 0
- cur_flowrate = df['flowrate'].tolist()
+ if nw_status == "off":
+ df.loc[:, "flowrate"] = 0
+ cur_flowrate = df["flowrate"].tolist()
return (True, True, Tout_val, cur_flowrate)
else:
# read Tout_val to dataframe
for i in range(n_BHE):
- df.loc[df.index[i], 'Tout_val'] = Tout_val[i]
+ df.loc[df.index[i], "Tout_val"] = Tout_val[i]
# TESPy solver
cur_Tin_val, cur_flowrate = get_tespy_results(t)
# check norm if network achieves the converge
if_success = False
pre_Tin_val = Tin_val
norm_dx = np.linalg.norm(
- abs(np.asarray(pre_Tin_val) - np.asarray(cur_Tin_val)))
+ abs(np.asarray(pre_Tin_val) - np.asarray(cur_Tin_val))
+ )
norm_x = np.linalg.norm(np.asarray(cur_Tin_val))
- if norm_dx/norm_x < 1e-6:
+ if norm_dx / norm_x < 1e-6:
if_success = True
# return to OGS
return (True, if_success, cur_Tin_val, cur_flowrate)
@@ -160,7 +171,7 @@ class BC(OpenGeoSys.BHENetwork):
# loading the TESPy model
project_dir = os.getcwd()
print("Project dir is: ", project_dir)
-nw = load_network('./pre/tespy_nw_closedloop')
+nw = load_network("./pre/tespy_nw_closedloop")
# set if print the network iteration info
nw.set_attr(iterinfo=False)
@@ -176,20 +187,20 @@ for i in range(n_BHE):
for c in nw.conns.index:
# bhe inlet and outlet conns
if c.target.label == data_index[i]: # inlet conns of bhe
- localVars['inlet_BHE' + str(i + 1)] = c
+ localVars["inlet_BHE" + str(i + 1)] = c
if c.source.label == data_index[i]: # outlet conns of bhe
- localVars['outlet_BHE' + str(i + 1)] = c
+ localVars["outlet_BHE" + str(i + 1)] = c
# time depended flowrate
-if switch_dyn_frate == 'on':
+if switch_dyn_frate == "on":
# import the name of inlet connection from the network csv file
- inlet_name = read_csv('./pre/tespy_nw/connections.csv',
- delimiter=';',
- index_col=[0]).iloc[0,0]
+ inlet_name = read_csv(
+ "./pre/tespy_nw/connections.csv", delimiter=";", index_col=[0]
+ ).iloc[0, 0]
for c in nw.conns.index:
# bhe inflow conns
if c.source.label == inlet_name: # inlet conns of bhe
- localVars['inlet_name'] = c
+ localVars["inlet_name"] = c
# instantiate BC objects referenced in OpenGeoSys
diff --git a/Tests/Data/Parabolic/T/3D_3BHEs_array/pre/3bhes.py b/Tests/Data/Parabolic/T/3D_3BHEs_array/pre/3bhes.py
index c5d32c0c14b..0d04ace81a6 100644
--- a/Tests/Data/Parabolic/T/3D_3BHEs_array/pre/3bhes.py
+++ b/Tests/Data/Parabolic/T/3D_3BHEs_array/pre/3bhes.py
@@ -7,8 +7,7 @@
# Execute this file to generate TESPy network csv files
from tespy.networks import network
-from tespy.components import (sink, source, splitter, merge,
- pump, heat_exchanger_simple)
+from tespy.components import sink, source, splitter, merge, pump, heat_exchanger_simple
from tespy.connections import connection, ref, bus
from tespy.tools.characteristics import char_line
from tespy.tools.data_containers import dc_cc
@@ -16,56 +15,59 @@ from tespy.tools.data_containers import dc_cc
import numpy as np
# %% network
-btes = network(fluids=['water'],
- T_unit='K',
- p_unit='bar',
- h_unit='kJ / kg',
- T_range=[273.25, 373.15],
- p_range=[1, 20],
- h_range=[1, 1000])
+btes = network(
+ fluids=["water"],
+ T_unit="K",
+ p_unit="bar",
+ h_unit="kJ / kg",
+ T_range=[273.25, 373.15],
+ p_range=[1, 20],
+ h_range=[1, 1000],
+)
# %% components
-fc_in = source('from consumer inflow')
-fc_out = sink('from consumer outflow')
+fc_in = source("from consumer inflow")
+fc_out = sink("from consumer outflow")
-pu = pump('pump')
+pu = pump("pump")
-sp = splitter('splitter', num_out=3)
+sp = splitter("splitter", num_out=3)
# bhe:
-bhe_name = 'BHE1'
-assert 'BHE1' in bhe_name, "BHE should be named with 'BHE1'"
+bhe_name = "BHE1"
+assert "BHE1" in bhe_name, "BHE should be named with 'BHE1'"
bhe1 = heat_exchanger_simple(bhe_name)
-bhe_name = 'BHE2'
-assert 'BHE2' in bhe_name, "BHE should be named with 'BHE2'"
+bhe_name = "BHE2"
+assert "BHE2" in bhe_name, "BHE should be named with 'BHE2'"
bhe2 = heat_exchanger_simple(bhe_name)
-bhe_name = 'BHE3'
-assert 'BHE3' in bhe_name, "BHE should be named with 'BHE3'"
+bhe_name = "BHE3"
+assert "BHE3" in bhe_name, "BHE should be named with 'BHE3'"
bhe3 = heat_exchanger_simple(bhe_name)
-mg = merge('merge', num_in=3)
+mg = merge("merge", num_in=3)
-cons = heat_exchanger_simple('consumer')
+cons = heat_exchanger_simple("consumer")
# %% connections
-fc_pu = connection(fc_in, 'out1', pu, 'in1')
+fc_pu = connection(fc_in, "out1", pu, "in1")
-pu_sp = connection(pu, 'out1', sp, 'in1')
+pu_sp = connection(pu, "out1", sp, "in1")
-sp_bhe1 = connection(sp, 'out1', bhe1, 'in1')
-sp_bhe2 = connection(sp, 'out2', bhe2, 'in1')
-sp_bhe3 = connection(sp, 'out3', bhe3, 'in1')
+sp_bhe1 = connection(sp, "out1", bhe1, "in1")
+sp_bhe2 = connection(sp, "out2", bhe2, "in1")
+sp_bhe3 = connection(sp, "out3", bhe3, "in1")
-bhe1_mg = connection(bhe1, 'out1', mg, 'in1')
-bhe2_mg = connection(bhe2, 'out1', mg, 'in2')
-bhe3_mg = connection(bhe3, 'out1', mg, 'in3')
+bhe1_mg = connection(bhe1, "out1", mg, "in1")
+bhe2_mg = connection(bhe2, "out1", mg, "in2")
+bhe3_mg = connection(bhe3, "out1", mg, "in3")
-mg_cons = connection(mg, 'out1', cons, 'in1')
+mg_cons = connection(mg, "out1", cons, "in1")
-cons_fc = connection(cons, 'out1', fc_out, 'in1')
+cons_fc = connection(cons, "out1", fc_out, "in1")
-btes.add_conns(fc_pu, pu_sp, sp_bhe1, sp_bhe2, sp_bhe3, bhe1_mg, bhe2_mg,
- bhe3_mg, mg_cons, cons_fc)
+btes.add_conns(
+ fc_pu, pu_sp, sp_bhe1, sp_bhe2, sp_bhe3, bhe1_mg, bhe2_mg, bhe3_mg, mg_cons, cons_fc
+)
# %% paramerization
@@ -73,25 +75,79 @@ btes.add_conns(fc_pu, pu_sp, sp_bhe1, sp_bhe2, sp_bhe3, bhe1_mg, bhe2_mg,
# pump
# flow_char
# provide volumetric flow in m^3 / s
-x = np.array([
- 0.00, 0.00001952885971862, 0.00390577194372, 0.005858657915586,
- 0.007811543887448, 0.00976442985931, 0.011717315831173, 0.013670201803035,
- 0.015623087774897, 0.017575973746759, 0.019528859718621, 0.021481745690483,
- 0.023434631662345, 0.025387517634207, 0.027340403606069, 0.029293289577931,
- 0.031246175549793, 0.033199061521655, 0.035151947493517, 0.037104833465379,
- 0.039057719437241, 0.041010605409104, 0.042963491380966, 0.044916377352828,
- 0.04686926332469, 0.048822149296552, 0.050775035268414, 0.052727921240276,
- 0.054680807212138, 0.056633693184
-])
+x = np.array(
+ [
+ 0.00,
+ 0.00001952885971862,
+ 0.00390577194372,
+ 0.005858657915586,
+ 0.007811543887448,
+ 0.00976442985931,
+ 0.011717315831173,
+ 0.013670201803035,
+ 0.015623087774897,
+ 0.017575973746759,
+ 0.019528859718621,
+ 0.021481745690483,
+ 0.023434631662345,
+ 0.025387517634207,
+ 0.027340403606069,
+ 0.029293289577931,
+ 0.031246175549793,
+ 0.033199061521655,
+ 0.035151947493517,
+ 0.037104833465379,
+ 0.039057719437241,
+ 0.041010605409104,
+ 0.042963491380966,
+ 0.044916377352828,
+ 0.04686926332469,
+ 0.048822149296552,
+ 0.050775035268414,
+ 0.052727921240276,
+ 0.054680807212138,
+ 0.056633693184,
+ ]
+)
# provide head in Pa
-y = np.array([
- 0.47782539, 0.47725723, 0.47555274, 0.47271192, 0.46873478, 0.46362130,
- 0.45737151, 0.44998538, 0.44146293, 0.43180416, 0.4220905, 0.40907762,
- 0.39600986, 0.38180578, 0.36646537, 0.34998863, 0.33237557, 0.31362618,
- 0.29374046, 0.27271841, 0.25056004, 0.22726535, 0.20283432, 0.17726697,
- 0.15056329, 0.12272329, 0.09374696, 0.06363430, 0.03238531, 0.00000000
-]) * 1e5
+y = (
+ np.array(
+ [
+ 0.47782539,
+ 0.47725723,
+ 0.47555274,
+ 0.47271192,
+ 0.46873478,
+ 0.46362130,
+ 0.45737151,
+ 0.44998538,
+ 0.44146293,
+ 0.43180416,
+ 0.4220905,
+ 0.40907762,
+ 0.39600986,
+ 0.38180578,
+ 0.36646537,
+ 0.34998863,
+ 0.33237557,
+ 0.31362618,
+ 0.29374046,
+ 0.27271841,
+ 0.25056004,
+ 0.22726535,
+ 0.20283432,
+ 0.17726697,
+ 0.15056329,
+ 0.12272329,
+ 0.09374696,
+ 0.06363430,
+ 0.03238531,
+ 0.00000000,
+ ]
+ )
+ * 1e5
+)
char = char_line(x=x, y=y)
pu.set_attr(flow_char=dc_cc(func=char, is_set=True))
@@ -105,8 +161,8 @@ bhe3.set_attr(D=0.013665, L=100, ks=0.00001)
# consumer
cons.set_attr(D=0.2, L=20, ks=0.00001)
# busses
-heat = bus('consumer heat demand')
-heat.add_comps({'c': cons, 'p': 'P'})
+heat = bus("consumer heat demand")
+heat.add_comps({"c": cons, "p": "P"})
btes.add_busses(heat)
# consumer heat demand
heat.set_attr(P=-3000) # W
@@ -116,7 +172,7 @@ heat.set_attr(P=-3000) # W
# system inlet
inflow_head = 2 # bar
-fc_pu.set_attr(p=inflow_head, m=0.6, fluid={'water': 1})
+fc_pu.set_attr(p=inflow_head, m=0.6, fluid={"water": 1})
# for BHEs:
# Tout:
@@ -128,8 +184,8 @@ bhe3_mg.set_attr(T=303.15)
pu_sp.set_attr(h=ref(cons_fc, 1, 0))
# %% solve
-btes.solve('design')
-#btes.print_results()
+btes.solve("design")
+# btes.print_results()
# %% save to csv:
-btes.save('tespy_nw', structure=True)
+btes.save("tespy_nw", structure=True)
diff --git a/Tests/Data/Parabolic/T/3D_3BHEs_array/pre/3bhes_closedloop.py b/Tests/Data/Parabolic/T/3D_3BHEs_array/pre/3bhes_closedloop.py
index 4ba9681b1cf..bce92ffb0cf 100644
--- a/Tests/Data/Parabolic/T/3D_3BHEs_array/pre/3bhes_closedloop.py
+++ b/Tests/Data/Parabolic/T/3D_3BHEs_array/pre/3bhes_closedloop.py
@@ -8,50 +8,112 @@
# Execute this file to generate TESPy network csv files
from tespy.networks import network
from tespy.connections import connection, ref
-from tespy.components import source, sink, pump, splitter, merge, heat_exchanger_simple, cycle_closer
+from tespy.components import (
+ source,
+ sink,
+ pump,
+ splitter,
+ merge,
+ heat_exchanger_simple,
+ cycle_closer,
+)
from tespy.tools import char_line, dc_cc
import numpy as np
# %% network
-btes = network(fluids=['water'], T_unit='K', p_unit='bar', h_unit='kJ / kg')
+btes = network(fluids=["water"], T_unit="K", p_unit="bar", h_unit="kJ / kg")
# %% components
-fc = cycle_closer('cycle closer')
-pu = pump('pump')
-sp = splitter('splitter', num_out=3)
+fc = cycle_closer("cycle closer")
+pu = pump("pump")
+sp = splitter("splitter", num_out=3)
# bhe:
-bhe1 = heat_exchanger_simple('BHE1')
-bhe2 = heat_exchanger_simple('BHE2')
-bhe3 = heat_exchanger_simple('BHE3')
+bhe1 = heat_exchanger_simple("BHE1")
+bhe2 = heat_exchanger_simple("BHE2")
+bhe3 = heat_exchanger_simple("BHE3")
-mg = merge('merge', num_in=3)
-cons = heat_exchanger_simple('consumer')
+mg = merge("merge", num_in=3)
+cons = heat_exchanger_simple("consumer")
## components paramerization
# pump
# flow_char
# provide volumetric flow in m^3 / s
-x = np.array([
- 0.00, 0.00001952885971862, 0.00390577194372, 0.005858657915586,
- 0.007811543887448, 0.00976442985931, 0.011717315831173, 0.013670201803035,
- 0.015623087774897, 0.017575973746759, 0.019528859718621, 0.021481745690483,
- 0.023434631662345, 0.025387517634207, 0.027340403606069, 0.029293289577931,
- 0.031246175549793, 0.033199061521655, 0.035151947493517, 0.037104833465379,
- 0.039057719437241, 0.041010605409104, 0.042963491380966, 0.044916377352828,
- 0.04686926332469, 0.048822149296552, 0.050775035268414, 0.052727921240276,
- 0.054680807212138, 0.056633693184
-])
+x = np.array(
+ [
+ 0.00,
+ 0.00001952885971862,
+ 0.00390577194372,
+ 0.005858657915586,
+ 0.007811543887448,
+ 0.00976442985931,
+ 0.011717315831173,
+ 0.013670201803035,
+ 0.015623087774897,
+ 0.017575973746759,
+ 0.019528859718621,
+ 0.021481745690483,
+ 0.023434631662345,
+ 0.025387517634207,
+ 0.027340403606069,
+ 0.029293289577931,
+ 0.031246175549793,
+ 0.033199061521655,
+ 0.035151947493517,
+ 0.037104833465379,
+ 0.039057719437241,
+ 0.041010605409104,
+ 0.042963491380966,
+ 0.044916377352828,
+ 0.04686926332469,
+ 0.048822149296552,
+ 0.050775035268414,
+ 0.052727921240276,
+ 0.054680807212138,
+ 0.056633693184,
+ ]
+)
# provide head in Pa
-y = np.array([
- 0.47782539, 0.47725723, 0.47555274, 0.47271192, 0.46873478, 0.46362130,
- 0.45737151, 0.44998538, 0.44146293, 0.43180416, 0.4220905, 0.40907762,
- 0.39600986, 0.38180578, 0.36646537, 0.34998863, 0.33237557, 0.31362618,
- 0.29374046, 0.27271841, 0.25056004, 0.22726535, 0.20283432, 0.17726697,
- 0.15056329, 0.12272329, 0.09374696, 0.06363430, 0.03238531, 0.00000000
-]) * 1e5
+y = (
+ np.array(
+ [
+ 0.47782539,
+ 0.47725723,
+ 0.47555274,
+ 0.47271192,
+ 0.46873478,
+ 0.46362130,
+ 0.45737151,
+ 0.44998538,
+ 0.44146293,
+ 0.43180416,
+ 0.4220905,
+ 0.40907762,
+ 0.39600986,
+ 0.38180578,
+ 0.36646537,
+ 0.34998863,
+ 0.33237557,
+ 0.31362618,
+ 0.29374046,
+ 0.27271841,
+ 0.25056004,
+ 0.22726535,
+ 0.20283432,
+ 0.17726697,
+ 0.15056329,
+ 0.12272329,
+ 0.09374696,
+ 0.06363430,
+ 0.03238531,
+ 0.00000000,
+ ]
+ )
+ * 1e5
+)
char = char_line(x=x, y=y)
pu.set_attr(flow_char=dc_cc(func=char, is_set=True))
pu.set_attr(eta_s=0.90)
@@ -64,29 +126,30 @@ bhe3.set_attr(D=0.013665, L=100, ks=0.00001)
# consumer
cons.set_attr(pr=1)
# consumer heat demand
-cons.set_attr(Q=-3000) # W
+cons.set_attr(Q=-3000) # W
# %% connections
-fc_pu = connection(fc, 'out1', pu, 'in1')
-pu_sp = connection(pu, 'out1', sp, 'in1')
+fc_pu = connection(fc, "out1", pu, "in1")
+pu_sp = connection(pu, "out1", sp, "in1")
-sp_bhe1 = connection(sp, 'out1', bhe1, 'in1')
-sp_bhe2 = connection(sp, 'out2', bhe2, 'in1')
-sp_bhe3 = connection(sp, 'out3', bhe3, 'in1')
+sp_bhe1 = connection(sp, "out1", bhe1, "in1")
+sp_bhe2 = connection(sp, "out2", bhe2, "in1")
+sp_bhe3 = connection(sp, "out3", bhe3, "in1")
-bhe1_mg = connection(bhe1, 'out1', mg, 'in1')
-bhe2_mg = connection(bhe2, 'out1', mg, 'in2')
-bhe3_mg = connection(bhe3, 'out1', mg, 'in3')
+bhe1_mg = connection(bhe1, "out1", mg, "in1")
+bhe2_mg = connection(bhe2, "out1", mg, "in2")
+bhe3_mg = connection(bhe3, "out1", mg, "in3")
-mg_cons = connection(mg, 'out1', cons, 'in1')
-cons_fc = connection(cons, 'out1', fc, 'in1')
+mg_cons = connection(mg, "out1", cons, "in1")
+cons_fc = connection(cons, "out1", fc, "in1")
-btes.add_conns(fc_pu, pu_sp, sp_bhe1, sp_bhe2, sp_bhe3, bhe1_mg, bhe2_mg,
- bhe3_mg, mg_cons, cons_fc)
+btes.add_conns(
+ fc_pu, pu_sp, sp_bhe1, sp_bhe2, sp_bhe3, bhe1_mg, bhe2_mg, bhe3_mg, mg_cons, cons_fc
+)
## connection parametrization
# system inlet
-fc_pu.set_attr(p=2, fluid={'water': 1})
+fc_pu.set_attr(p=2, fluid={"water": 1})
# for BHEs:
# Tout:
@@ -96,8 +159,8 @@ bhe3_mg.set_attr(T=303.15)
# %% solve
-btes.solve('design')
-#btes.print_results()
+btes.solve("design")
+# btes.print_results()
# %% save to csv:
-btes.save('tespy_nw_closedloop')
+btes.save("tespy_nw_closedloop")
diff --git a/Tests/Data/ThermoMechanics/BDT/generate_ref.py b/Tests/Data/ThermoMechanics/BDT/generate_ref.py
index 9b694cfd2d0..b94584435dc 100755
--- a/Tests/Data/ThermoMechanics/BDT/generate_ref.py
+++ b/Tests/Data/ThermoMechanics/BDT/generate_ref.py
@@ -12,9 +12,9 @@ grid = reader.GetOutput()
nnodes = grid.GetNumberOfPoints()
ref_A = np.loadtxt("bdt.res")
-ts = ref_A[:,0]
-epss = ref_A[:,1:(1+6)]
-sigmas = ref_A[:,(1+6):(1+12)]
+ts = ref_A[:, 0]
+epss = ref_A[:, 1 : (1 + 6)]
+sigmas = ref_A[:, (1 + 6) : (1 + 12)]
for t, eps, sigma in zip(ts, epss, sigmas):
# assert that we don't have to bother with permuting the shear
--
GitLab