Merge pull request #54 from mattfalcone1997/nonconformal-tests

Author: mattfalcone1997

.gitignore

@@ -145,6 +145,12 @@ moose/include/base/Precompiled.h.gch
 .*.swp
 .*.swo
 
+# clang
+.clangd
+
+# OpenFOAM mesh object files
+test/**/*.obj
+
 framework/contrib/exodiff/exodiff
 
 # Mac garbage

run_tests

@@ -8,5 +8,11 @@ MOOSE_DIR = os.environ.get('MOOSE_DIR', MOOSE_DIR)
 
 sys.path.append(os.path.join(MOOSE_DIR, 'python'))
 
+# add python test to PYTHONPATH
+HIPPO_TEST_PYTHON_DIR = os.path.abspath(os.path.join(os.path.dirname(__file__), 'test', 'python'))
+pythonpath = os.environ.get('PYTHONPATH', '').split(':')
+pythonpath = [HIPPO_TEST_PYTHON_DIR] + pythonpath
+os.environ['PYTHONPATH'] = ':'.join(pythonpath)
+
 from TestHarness import TestHarness
 TestHarness.buildAndRun(sys.argv, 'hippo', MOOSE_DIR)

test/python/analytical.py

@@ -0,0 +1,97 @@
+""" # analytical.py
+
+Module contains analytical solutions used in the test scripts
+"""
+
+import numpy as np
+from scipy.special import erfc
+
+def unsteady1d_temp(x: np.ndarray,
+                    time: float,
+                    temp_cold: float,
+                    temp_hot: float,
+                    k1: float,
+                    k2: float,
+                    rho_cp1: float,
+                    rho_cp2: float) -> np.ndarray:
+    """Returns the analytical solution to 1D unsteady heat conduction problem.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        array of x coordinates to evaluate solution
+    time : float
+        time to evaluate solution
+    temp_cold : float
+        Cold temperature
+    temp_hot : float
+        Hot temperautre
+    k1 : float
+        thermal conductivity for x<0
+    k2 : float
+        thermal conductivity for x>0
+    rho_cp1 : float
+        \\rho c_p for  x<0
+    rho_cp2 : float
+        \\rho c_p for  x>0
+
+    Returns
+    -------
+    np.ndarray
+        Temperature for x at time
+    """
+
+    def temp_minus(
+        x: np.ndarray,
+        time: float,
+        temp_cold: float,
+        temp_hot: float,
+        k1: float,
+        k2: float,
+        rho_cp1: float,
+        rho_cp2: float,
+    ) -> np.ndarray:
+        numerator = (temp_hot - temp_cold) * np.sqrt(k2 * rho_cp2)
+        denominator = np.sqrt(k1 * rho_cp1) + np.sqrt(k2 * rho_cp2)
+        erfc_arg = -x / (2 * np.sqrt(k1 / rho_cp1 * time))
+        return temp_hot - (numerator / denominator) * erfc(erfc_arg)
+
+    def temp_plus(
+        x: np.ndarray,
+        time: float,
+        temp_cold: float,
+        temp_hot: float,
+        k1: float,
+        k2: float,
+        rho_cp1: float,
+        rho_cp2: float,
+    ) -> np.ndarray:
+        numerator = (temp_hot - temp_cold) * np.sqrt(k1 * rho_cp1)
+        denominator = np.sqrt(k1 * rho_cp1) + np.sqrt(k2 * rho_cp2)
+        erfc_arg = x / (2 * np.sqrt(k2 / rho_cp2 * time))
+        return temp_cold + (numerator / denominator) * erfc(erfc_arg)
+
+    positive_x_idx = x >= 0
+    temp = np.empty_like(x)
+    temp[positive_x_idx] = temp_plus(
+        x=x[positive_x_idx],
+        time=time,
+        temp_cold=temp_cold,
+        temp_hot=temp_hot,
+        k1=k1,
+        k2=k2,
+        rho_cp1=rho_cp1,
+        rho_cp2=rho_cp2,
+    )
+    negative_x_idx = ~positive_x_idx
+    temp[negative_x_idx] = temp_minus(
+        x=x[negative_x_idx],
+        time=time,
+        temp_cold=temp_cold,
+        temp_hot=temp_hot,
+        k1=k1,
+        k2=k2,
+        rho_cp1=rho_cp1,
+        rho_cp2=rho_cp2,
+    )
+    return temp

test/python/read_hippo_data.py

@@ -0,0 +1,54 @@
+
+from pathlib import Path
+from typing import Literal
+import numpy as np
+import pyvista as pv
+
+def read_moose_exodus_data(exo_file,
+                        time,
+                        variable,
+                        block: int | str=0):
+
+    reader: pv.ExodusIIReader = pv.get_reader(exo_file)
+    reader.set_active_time_value(time)
+
+    if block != 0 or block != 'Element blocks':
+        reader.node_sets.enable_all()
+        reader.side_sets.enable_all()
+        reader.node_sets.enable_all_arrays()
+        reader.side_sets.enable_all_arrays()
+
+    data: pv.UnstructuredGrid = reader.read()[block].combine(True)
+
+    if variable in data.point_data.keys():
+        coords = data.points
+        foam_variable = data.point_data[variable]
+    elif variable in data.cell_data.keys():
+        coords = data.cell_centers().points
+        foam_variable = data.cell_data[variable]
+
+    else:
+        raise KeyError(f"Variable {variable} not found")
+
+    foam_coords = dict(zip(('x', 'y', 'z'), coords.T))
+    return foam_coords, foam_variable
+
+def read_openfoam_data(case_dir: Path | str | bytes,
+                       time: float,
+                       variable: str,
+                       block: str = 'internalMesh',
+                       case_type: Literal['decomposed', 'reconstructed'] = 'reconstructed'
+                       ) -> tuple[np.ndarray, np.ndarray]:
+
+    file_name = Path(case_dir) / "case.foam"
+    file_name.touch()
+
+    reader: pv.POpenFOAMReader = pv.get_reader(file_name)
+    reader.case_type = case_type
+    reader.set_active_time_value(time)
+    data: pv.UnstructuredGrid = reader.read()[block]
+
+    coords = data.cell_centers().points
+    foam_variable = data.cell_data[variable]
+    foam_coords = dict(zip(('x', 'y', 'z'), coords.T))
+    return foam_coords, foam_variable

test/tests/mesh/polygonal/gold/run_temp_over_line_csv_temp_over_line_0032.csv

@@ -42,7 +42,7 @@
     []
 
     [heat_flux_to_fluid]
-        type = MultiAppProjectionTransfer
+        type = MultiAppGeneralFieldShapeEvaluationTransfer
         source_variable = wall_heat_flux
         to_multi_app = hippo
         variable = solid_heat_flux

test/tests/mesh/polygonal/run.i

@@ -0,0 +1,51 @@
+"""Test module for polygonal OpenFOAM meshes"""
+
+import unittest
+from pathlib import Path
+
+import numpy as np
+
+#import Hippo test python functions
+from analytical import unsteady1d_temp
+from read_hippo_data import read_moose_exodus_data, read_openfoam_data
+
+
+RUN_DIR = Path(__file__).parent
+FOAM_CASE = RUN_DIR / "fluid-openfoam"
+K_SOLID = 1
+RHO_CP_SOLID = 1
+K_FLUID = 4
+RHO_CP_FLUID = 16
+T_HOT = 2.
+T_COLD = 1.
+L = 1
+
+
+class TestUnsteadyHeatConductionInInfiniteSystem(unittest.TestCase):
+    """Test class for 1D unsteady problem with polygonal elements for OpenFOAM"""
+    def test_matches_analytic_solution_at_times(self):
+        """Compare against analytical solution."""
+        times = [0.0025, 0.005, 0.01]  # seconds
+        for time in times:
+            moose_coords, moose_temperature = read_moose_exodus_data(RUN_DIR / "run_out.e",
+                                                                time,
+                                                                "temp")
+            foam_coords, foam_temperature = read_openfoam_data(FOAM_CASE,
+                                                          time,
+                                                          'T')
+            x = np.concatenate([moose_coords['x'], foam_coords['x']])
+            temp = np.concatenate([moose_temperature, foam_temperature])
+
+            analytic_temp = unsteady1d_temp(
+                x=x,
+                time=time,
+                temp_cold=T_COLD,
+                temp_hot=T_HOT,
+                k1=K_SOLID,
+                k2=K_FLUID,
+                rho_cp1=RHO_CP_SOLID,
+                rho_cp2=RHO_CP_FLUID,
+            )
+
+            rmse = np.sqrt(np.sum(np.square(analytic_temp - temp)) / len(temp))
+            self.assertLess(rmse, 5e-3, msg=f"for time = {time} s")

test/tests/mesh/polygonal/test.py

@@ -2,13 +2,18 @@
   [poly]
     [setup]
       type = RunCommand
-      command = 'bash -c "cd fluid-openfoam && ./Allclean && fluentMeshToFoam constant/polygon.msh && polyDualMesh -overwrite 35"'
+      command = 'bash -c "cd fluid-openfoam && ./Allclean && fluentMeshToFoam constant/polygon.msh && polyDualMesh -overwrite 75"'
     []
-    [verify]
-      type = CSVDiff
+    [run]
+      type = RunApp
       input = run.i
-      csvdiff = run_temp_over_line_csv_temp_over_line_0032.csv
       prereq = poly/setup
     []
+
+    [verify]
+      type = PythonUnitTest
+      input = test.py
+      prereq = poly/run
+    []
   []
 []