Hi,
I am trying to adapt the openfoam tutorial flow-over-heated-plate to use a in house solver called TMNGSolve to simmulate the solid component.
I get almost identical simulation results to the openFOAM-openFOAM refrence case, except for two differences:
-
The simulations disagree in the first two timesteps
-
From the third timestep, the TMNGSolve simmulation seems to be one timestep behind.
To illustrate this fact i plotted the average temperature on the interface for the first 10 timesteps in the simmulation. As you can se the average from the openfoam simmulation for timestep 5 equals the average in the TMNGSolve simmulation for timestep 6.
Can you see any mistake in the solid-TMNGSolve.py file wich can explain this difference?
(i include the plot of the average surface temperature, precice-config, solidTMNGSolve.py and the code i used to make the plot)
Best regards,
Eirik Høydalsvik
Blockquote
solid-TMNGSOLVE:
from TMNGSolve_setup import update_boundary_and_solve, solver_args, top_nodes, top_pos, get_conv_arr, T_bottom, Flux_zero, l, w
import numpy as np
import precice
import matplotlib.pyplot as plt
import pandas as pd
import os
import sys
import ngsolve as ng
import copy
path2add = os.path.normpath( # nopep8
os.path.abspath('/home/eirikhoy/code_dir/tmngsolve/tmngsolve/')) # nopep8
if (not (path2add in sys.path)): # nopep8
sys.path.append(path2add) # nopep8
from utils import get_fluxes_projection # nopep8
from solvers import solve_transient_heat_equation # nopep8
# *************
# setup preCICE
x = top_pos
nx = len(x)
T0 = 310 # Kelvin init temp
# Setter opp array for å ta imot temperatur fra OpenFOAM. Første omgang bare 1D array.
T = np.ones(nx)*T0
participant_name = "Solid"
# Definerer hvilken participant dette er, samt hvordan fil de skal "kommunisere" gjennom
config_file_name = "../precice-config.xml"
solver_process_index = 0
solver_process_size = 1
interface = precice.Interface(
participant_name, config_file_name, solver_process_index, solver_process_size)
mesh_name = "Solid-Mesh" # Definerer meshet.
mesh_id = interface.get_mesh_id(mesh_name)
temp_name = "Temperature"
# De parametrene som skal brukes, se precice-config.xml hvordan de kommuniserer.
heatFlux_name = "Heat-Flux"
temp_id = interface.get_data_id(temp_name, mesh_id)
heatFlux_id = interface.get_data_id(heatFlux_name, mesh_id)
vertices = [[x0, 0] for x0 in x]
vertex_ids = interface.set_mesh_vertices(mesh_id, vertices)
dt = 0.01
T_top_init_arr = np.array(
[top_nodes, solver_args["ic"] * np.ones_like(top_nodes)]).T
precice_dt = interface.initialize()
t = 0
num_sol = solver_args["ic"]
# setting up variables for writing heat flux
q = np.zeros(nx)
y = l*np.ones(nx) # w is length of rectangle in y-dir
while interface.is_coupling_ongoing():
# write checkpoint
if interface.is_action_required(precice.action_write_iteration_checkpoint()):
num_sol_checkpoint = num_sol
t_checkpoint = t
interface.mark_action_fulfilled(
precice.action_write_iteration_checkpoint())
solver_args["ic"] = num_sol
dt = np.minimum(dt, precice_dt)
T = interface.read_block_scalar_data(
temp_id, vertex_ids) # leser temperatur inn i T-array
# Updating boundary conditions and time step in solver_args
num_sol, q = update_boundary_and_solve(solver_args, x, y, dt, T, num_sol)
precice_dt = interface.advance(dt) # Sender/mottar, går til neste tid
t += dt
interface.write_block_scalar_data(heatFlux_id, vertex_ids, q)
# read checkpoint
if interface.is_action_required(precice.action_read_iteration_checkpoint()):
num_sol = num_sol_checkpoint
t = t_checkpoint
interface.mark_action_fulfilled(
precice.action_read_iteration_checkpoint())
interface.finalize()
vtk = ng.VTKOutput(ma=solver_args["mesh"],
coefs=[num_sol],
names=["Temperature"],
filename=f"output/final",
subdivision=0)
vtk.Do()
Precice config:
<?xml version="1.0"?>
<precice-configuration>
<log>
<sink
filter="%Severity% > debug and %Rank% = 0"
format="---[precice] %ColorizedSeverity% %Message%"
enabled="true" />
</log>
<solver-interface dimensions="2">
<data:scalar name="Temperature" />
<data:scalar name="Heat-Flux" />
<mesh name="Fluid-Mesh">
<use-data name="Temperature" />
<use-data name="Heat-Flux" />
</mesh>
<mesh name="Solid-Mesh">
<use-data name="Temperature" />
<use-data name="Heat-Flux" />
</mesh>
<participant name="Fluid">
<use-mesh name="Fluid-Mesh" provide="yes" />
<use-mesh name="Solid-Mesh" from="Solid" />
<read-data name="Heat-Flux" mesh="Fluid-Mesh" />
<write-data name="Temperature" mesh="Fluid-Mesh" />
<mapping:nearest-neighbor
direction="read"
from="Solid-Mesh"
to="Fluid-Mesh"
constraint="consistent" />
</participant>
<participant name="Solid">
<export:vtk directory="preCICE-output" />
<use-mesh name="Fluid-Mesh" from="Fluid" />
<use-mesh name="Solid-Mesh" provide="yes" />
<mapping:nearest-neighbor
direction="read"
from="Fluid-Mesh"
to="Solid-Mesh"
constraint="consistent" />
<read-data name="Temperature" mesh="Solid-Mesh" />
<write-data name="Heat-Flux" mesh="Solid-Mesh" />
</participant>
<m2n:sockets from="Fluid" to="Solid" exchange-directory=".." />
<coupling-scheme:serial-implicit>
<time-window-size value="0.01" />
<max-time value="1" />
<max-iterations value="200" />
<participants first="Fluid" second="Solid" />
<exchange data="Temperature" mesh="Fluid-Mesh" from="Fluid" to="Solid" />
<exchange data="Heat-Flux" mesh="Solid-Mesh" from="Solid" to="Fluid" />
<relative-convergence-measure limit="1.0e-5" data="Temperature" mesh="Fluid-Mesh" />
<acceleration:aitken>
<data mesh="Solid-Mesh" name="Heat-Flux" />
<initial-relaxation value="0.5" />
</acceleration:aitken>
</coupling-scheme:serial-implicit>
</solver-interface>
</precice-configuration>
Plotting code:
import vtk
from matplotlib import pyplot as plt
import numpy as np
def vtk_to_dict(case, time):
vtkFileName = "solid-{}/preCICE-output/Fluid-Mesh-Solid.{}.vtk".format(
case, time)
# read the vtk file as an unstructured grid
reader = vtk.vtkUnstructuredGridReader()
reader.SetFileName(vtkFileName)
reader.ReadAllVectorsOn()
reader.ReadAllScalarsOn()
reader.Update()
# obtain the data
data = reader.GetOutput()
n_data = data.GetPointData().GetNumberOfTuples()
name = "Temperature"
data_names = []
i = 0
data_id = 0
max_i = data.GetPointData().GetNumberOfArrays()
while i < max_i:
this_data_name = data.GetPointData().GetArray(i).GetName()
data_names.append(this_data_name)
if (this_data_name == name):
data_id = i
break
i += 1
data_dict = {}
if not data_id:
raise Exception(
"For file {} name {} not found. Only the following names are available: {}. "
"Aborting!".format(vtkFileName, name, data_names))
for i in range(n_data):
data_dict[data.GetPoint(i)] = data.GetPointData().GetArray(
data_id).GetValue(i)
return data_dict
cases = []
# cases.append('python')
cases.append('openfoam')
cases.append('TMNGSolve')
case_labels = {'fenics': 'OpenFOAM-FEniCS', 'openfoam': 'OpenFOAM-OpenFOAM',
'nutils': 'OpenFOAM-Nutils', 'TMNGSolve': 'OpenFOAM-TMNGSolve', 'python': 'OpenFOAM-HTC'}
styles = [':', '-', '--']
colors = ['r', 'b', 'g', 'k']
def plot_surface(time1, time2, savefile):
fig, ax = plt.subplots()
data_dict = {}
for i, case in enumerate(cases):
case_data1 = vtk_to_dict(case, time1)
x1, t1 = [p[0] for p in case_data1.keys()], np.array(
list(case_data1.values()))
# theta = (t1 - 300) / (310 - 300)
data_dict[case] = {"x": np.array(x1)}
data_dict[case]["t"] = np.array(t1)
ax.plot(np.array(x1), t1, colors[i % 4] + styles[i %
3], label=f"{case_labels[case]}{time1}")
case_data2 = vtk_to_dict(case, time2)
x2, t2 = [p[0] for p in case_data2.keys()], np.array(
list(case_data2.values()))
plt.plot(x2, t2, colors[(i+2) % 4] + styles[i %
3], label=f"{case_labels[case]}{time2}")
print(data_dict['openfoam']['x'] - data_dict['TMNGSolve']['x'])
delT = data_dict['openfoam']['t'] - data_dict['TMNGSolve']['t']
Topenfoam = data_dict['openfoam']['t']
print(np.mean(np.abs(delT/Topenfoam)))
FONTSIZE = 14
ax.set_ylabel("T (K)", fontsize=FONTSIZE)
ax.set_xlabel("x (m)", fontsize=FONTSIZE)
ax.legend(fontsize=FONTSIZE)
plt.savefig(f"plots/{savefile}.pdf")
plt.tight_layout()
plt.show()
def plot_average_interface(t1, t2, savefile):
time_l = np.linspace(t1, t2, t2-t1+1)
str_l = [f"dt{int(t)}" for t in time_l]
surface_average_dict = {case[1]: np.zeros_like(
time_l) for case in enumerate(cases)}
for i, time_str in enumerate(str_l):
for j, case in enumerate(cases):
case_data = vtk_to_dict(case, time_str)
x1, t = [p[0] for p in case_data.keys()], np.array(
list(case_data.values()))
print(case, np.average(t))
surface_average_dict[case][i] = np.average(t)
fig, ax = plt.subplots()
for case in surface_average_dict:
ax.plot(time_l, surface_average_dict[case], ".", label=case)
ax.legend()
ax.set_xlabel("timestep")
ax.set_ylabel("average surface temperature (C)")
ax.set_xticks(time_l)
plt.savefig(f"plots/{savefile}.pdf")
plt.show()
plot_average_interface(1, 10, "average_through_time")
