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
average_through_time

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")
1 Like

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:

  1. read T
  2. solve
  3. advance
  4. write q

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

  1. read T
  2. solve
  3. write q
  4. advance

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

How does this change the results?

Thank you for your response:-)

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”:

average_through_time_serial_implicit_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:

average_through_time_serial_implicit_Fluid-Mesh-Solid

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:
average_through_time_serial_implicit_Fluid-Mesh-Solid

average_through_time_serial_implicit_Solid-Mesh-Solid

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:

average_through_time_parallel_implicit_Fluid-Mesh-Solid
average_through_time_parallel_implicit_Solid-Mesh-Solid

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