Solid solver implemented with python API and in-house solver is one time-step behind openFOAM-openFOAM refrence case

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:

1. The simulations disagree in the first two timesteps

2. 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

os.path.abspath('/home/eirikhoy/code_dir/tmngsolve/tmngsolve/'))  # nopep8
if (not (path2add in sys.path)):  # 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)

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)

num_sol = num_sol_checkpoint
t = t_checkpoint

interface.mark_action_fulfilled(

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" />
<write-data name="Temperature" mesh="Fluid-Mesh" />
<mapping:nearest-neighbor
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
from="Fluid-Mesh"
to="Solid-Mesh"
constraint="consistent" />
<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

# obtain the data
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")
``````

I would not be surprised if the time stepping in the OpenFOAM adapter has issues, but I would start first by reconsidering the order of events in this code:

it essentially does:

2. solve
4. write q

However, the `advance()` is exchanging old information. It should, instead, be:

2. solve
3. write q

So, `advance` needs to have updated buffers to send updated values to the other participant.

How does this change the results?

I simply changed the order of the advance and write q like you said:

``````    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)

# skrive varmefluks til preCICE
interface.write_block_scalar_data(heatFlux_id, vertex_ids, q)

precice_dt = interface.advance(dt)  # Sender/mottar, går til neste tid
t += dt
``````

And got the following results when plotting from “Solid-Mesh-Solid”:

The results are a lot better so your suggestion seems to help:-)

However, when i plot from “Fluid-Mesh-Solid” i get almost identical results from the two simmulations become almost identical:

Can you explain the dfifference between Solid-Mesh-Solid and Fluid-Mesh-Solid?

Can you also explain what you mean by “updated buffers”? Naively i would think that the order of advance and write would not matter since write does not use the “precice_dt” variable.

I assume you are using a serial-implicit scheme. It would be interesting to try a parallel-implicit scheme, or make the convergence criteria tighter.

The second participant is probably behind the first one, because in an implicit scheme, the coupling iterates, updating the values computed.

The `precice_dt` is not relevant here. When calling `advance`, preCICE exchanges the data that has been passed to `write_block_scalar_data`, which needs to be updated first. By “buffer”, I mean here the preCICE data that the `q` variable is associated with. You update that data using `write_block_scalar_data`.

I changed the relative convergence measure from 1e-5 to 1e-7 and it did not change the results.

``````<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-7" 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>
``````

Plots:

When it comes to paralell-implicit, i used the following configuration:

``````<coupling-scheme:parallel-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-7" data="Temperature" mesh="Fluid-Mesh" />
<acceleration:IQN-ILS>
<data mesh="Solid-Mesh" name="Heat-Flux" />
<preconditioner type="residual-sum"/>
<filter type="QR2" limit="1e-3"/>
<initial-relaxation value="0.1"/>
<max-used-iterations value="100"/>
<time-windows-reused value="20"/>
</acceleration:IQN-ILS>
</coupling-scheme:parallel-implicit>
``````

But then, none of the simmulations gave the same results:

Sorry for the delay, @EirikJaccheri, too much going on lately…

Did you manage to understand/fix this issue, or should I try to look deeper?

Yes, i have fixed the issue. Thanks for your help:-)