diff --git a/apps/tutorials/codegen/01_CodegenHeatEquation.cpp b/apps/tutorials/codegen/01_CodegenHeatEquation.cpp
index 361ebf76434fa47388b6fe896f7f9458c3b72929..8f425ca008f8b4609977b50573943ef5d2c53d56 100644
--- a/apps/tutorials/codegen/01_CodegenHeatEquation.cpp
+++ b/apps/tutorials/codegen/01_CodegenHeatEquation.cpp
@@ -67,7 +67,7 @@ void initDirichletBoundaryNorth(shared_ptr< StructuredBlockForest > blocks, Bloc
{
const Vector3< real_t > p = blocks->getBlockLocalCellCenter(*block, *cell);
// Set the dirichlet boundary to f(x) = 1 + sin(x) * x^2
- real_t v = real_c(1.0 + std::sin(2 * math::pi * p[0]) * p[0] * p[0]);
+ real_t v = real_c(2000.0 + 1500 * std::sin(2 * math::pi * p[0]) * p[0] * p[0]);
u->get(*cell) = v;
u_tmp->get(*cell) = v;
}
@@ -103,7 +103,7 @@ int main(int argc, char** argv)
WALBERLA_CHECK_FLOAT_EQUAL(dx, dy);
- const real_t dt = real_c(1e-4);
+ const real_t dt = real_c(1);
const real_t kappa = real_c(1.0);
///////////////////////////
@@ -120,8 +120,9 @@ int main(int argc, char** argv)
/// FIELDS ///
//////////////
- BlockDataID uFieldId = field::addToStorage< ScalarField >(blocks, "u", real_c(0.0), field::fzyx, uint_c(1));
- BlockDataID uTmpFieldId = field::addToStorage< ScalarField >(blocks, "u_tmp", real_c(0.0), field::fzyx, uint_c(1));
+ BlockDataID uFieldId = field::addToStorage< ScalarField >(blocks, "u", real_c(2000.0), field::fzyx, uint_c(1));
+ BlockDataID uTmpFieldId = field::addToStorage< ScalarField >(blocks, "u_tmp", real_c(2000.0), field::fzyx, uint_c(1));
+ BlockDataID kappa_outFieldId = field::addToStorage< ScalarField >(blocks, "kappa_out", real_c(0.0), field::fzyx, uint_c(1));
/////////////////////
/// COMMUNICATION ///
@@ -153,13 +154,16 @@ int main(int argc, char** argv)
SweepTimeloop timeloop(blocks, uint_c(2e4));
timeloop.add() << BeforeFunction(commScheme, "Communication") << BeforeFunction(neumann, "Neumann Boundaries")
- << Sweep(pystencils::HeatEquationKernel(uFieldId, uTmpFieldId, dt, dx, kappa), "HeatEquationKernel")
+ << Sweep(pystencils::HeatEquationKernel(kappa_outFieldId, uFieldId, uTmpFieldId, dt, dx, kappa), "HeatEquationKernel")
<< AfterFunction([blocks, uFieldId, uTmpFieldId]() { swapFields(*blocks, uFieldId, uTmpFieldId); },
"Swap");
auto vtkWriter = field::createVTKOutput< ScalarField, float >( uFieldId, *blocks, "temperature", uint_c(200), uint_c(0) );
vtkWriter();
timeloop.addFuncAfterTimeStep(vtkWriter, "VTK");
+ auto vtkWriter1 = field::createVTKOutput< ScalarField, float >( kappa_outFieldId, *blocks, "kappa", uint_c(200), uint_c(0) );
+ vtkWriter1();
+ timeloop.addFuncAfterTimeStep(vtkWriter1, "VTK1");
timeloop.run();
diff --git a/apps/tutorials/codegen/HeatEquationKernel.py b/apps/tutorials/codegen/HeatEquationKernel.py
index 024940a38dc1f75e748db4c1987800a9a85cb178..9c418fc977b34f24344619c04e40d565a2210980 100644
--- a/apps/tutorials/codegen/HeatEquationKernel.py
+++ b/apps/tutorials/codegen/HeatEquationKernel.py
@@ -1,27 +1,278 @@
+import sys
+sys.path.append('/local/ca00xebo/repos/walberla')
+import numpy as np
import sympy as sp
import pystencils as ps
from pystencils_walberla import CodeGeneration, generate_sweep
+from dataclasses import dataclass
+from pystencils.typing import PointerType
+import src.materials.alloys.Ti6Al4V as ti
+from src.materials.get_parameters import get_parameters
+from src.materials.alloys.Ti6Al4V import Ti6Al4V_parameters
+'''
+def interpolate(x, xmin, xmax, ymin, ymax):
+ """
+ Linear interpolation between two points: y = y_1 + ((y2 - y1) / (x2 - x1)) * (x - x1)
+ """
+ return (ymin * (xmax - x) + ymax * (x - xmin)) / (xmax - xmin)
+
+
+def interpolate_property(T, base_temperature, base_values): # TODO
+ """
+ Interpolates properties based on temperature for both symbolic and numeric cases.
+
+ :param T: Temperature value, either symbolic or numeric.
+ :param base_temperature: List or array of base temperatures.
+ :param base_values: List or array of property values corresponding to base temperatures.
+ """
+ # Ensure base_temperature and base_values are lists or arrays
+ if not isinstance(base_temperature, (list, np.ndarray)) or not isinstance(base_values, (list, np.ndarray)):
+ raise TypeError("base_temperature and base_values must be of type list or np.ndarray")
+ # Convert to lists if they are numpy arrays
+ if isinstance(base_temperature, np.ndarray):
+ base_temperature = base_temperature.tolist()
+ if isinstance(base_values, np.ndarray):
+ base_values = base_values.tolist()
+ # Check if base_temperature and base_values have the same length
+ if len(base_temperature) != len(base_values):
+ raise ValueError("base_temperature and base_values must have the same length")
+ # assert len(base_temperature) == len(base_values), "base_temperature and base_values must have the same length"
+
+ if isinstance(T, sp.Expr): # Symbolic expression
+ base_temperature = [sp.Float(t) for t in base_temperature]
+ base_values = [sp.Float(val) for val in base_values]
+ conditions = [
+ (base_values[0], T < base_temperature[0]),
+ (base_values[-1], T >= base_temperature[-1])]
+ for i in range(len(base_temperature) - 1):
+ interp_expr = (base_values[i] +
+ (base_values[i+1] - base_values[i]) / (base_temperature[i+1] - base_temperature[i]) *
+ (T - base_temperature[i]))
+ conditions.append((interp_expr, sp.And(T >= base_temperature[i], T < base_temperature[i + 1])))
+ return sp.Piecewise(*conditions)
+ else: # Numeric expression
+ return np.interp(T, base_temperature, base_values)
+
+
+def density_by_thermal_expansion(T, base_temperature, base_density, thermal_expansion_coefficient):
+ """
+ Calculate density based on thermal expansion: rho_0 / (1 + tec * (T - T_0))
+
+ :param T: Temperature value.
+ :param base_temperature: Reference temperature.
+ :param base_density: Density at reference temperature.
+ :param thermal_expansion_coefficient: Thermal expansion coefficient.
+ """
+ return base_density / (1 + thermal_expansion_coefficient * (T - base_temperature)) ** 3
+
+
+def thermal_diffusivity_func(heat_conductivity, density, specific_heat_capacity):
+ """
+ Calculate thermal diffusivity: k(T) / (rho(T) * c_p(T))
+
+ :param heat_conductivity: Thermal conductivity.
+ :param density: Density.
+ :param specific_heat_capacity: Specific heat capacity.
+ """
+ return heat_conductivity / (density * specific_heat_capacity)
+
+
+def density_lookup(T, temperature_array, density_array):
+ """
+ Look up density based on temperature using interpolation.
+
+ :param T: Temperature value.
+ :param temperature_array: Array of temperatures.
+ :param density_array: Array of densities.
+ """
+ if isinstance(T, sp.Expr): # Symbolic expression
+ conditions = [
+ (density_array[0], T < temperature_array[0]),
+ (density_array[-1], T >= temperature_array[-1])
+ ]
+ for i in range(len(temperature_array) - 1):
+ interp_expr = (density_array[i] +
+ (density_array[i+1] - density_array[i]) / (temperature_array[i+1] - temperature_array[i]) *
+ (T - temperature_array[i]))
+ conditions.append((interp_expr, sp.And(T >= temperature_array[i], T < temperature_array[i + 1])))
+ return sp.Piecewise(*conditions)
+ else: # Numeric expression
+ return np.interp(T, temperature_array, density_array)
+
+
+def get_arguments(T, temperatures, args):
+ """
+ Recursively get arguments for parameter functions.
+
+ :param T: Temperature value.
+ :param temperatures: Solidus and liquidus temperatures.
+ :param args: Arguments to be processed.
+ """
+ converted_args = []
+ for arg in args:
+ if isinstance(arg, (Solidus_Liquidus_Parameter, Functional_Parameter)):
+ converted_args.append(get_parameters(T, temperatures, arg))
+ else:
+ converted_args.append(arg)
+ return converted_args
+
+
+def get_parameters(T, temperatures, parameter):
+ """
+ Resolve parameters based on temperature and parameter type.
+
+ :param T: Temperature value.
+ :param temperatures: Solidus and liquidus temperatures.
+ :param parameter: Parameter to be resolved.
+ """
+ if isinstance(parameter, Solidus_Liquidus_Parameter):
+ return sp.Piecewise(
+ (parameter.sol, T <= temperatures.sol),
+ (parameter.liq, T >= temperatures.liq),
+ (interpolate(T, temperatures.sol, temperatures.liq, parameter.sol, parameter.liq), True)
+ )
+ elif isinstance(parameter, Functional_Parameter):
+ resolved_args = [get_parameters(T, temperatures, arg) for arg in parameter.args]
+ if parameter.func == 'thermal_expansion':
+ return density_by_thermal_expansion(T, *resolved_args) # *get_parameters(T, temperatures, parameter.args)
+ elif parameter.func == 'interpolation':
+ return interpolate_property(T, *resolved_args)
+ elif parameter.func == 'interpolation_in_array': # TODO
+ T_base, T_incr, data_array, data_symbol = resolved_args
+ idx = sp.floor((T - T_base) / T_incr)
+ mod = (T - T_base) - idx * T_incr # mod replaces (T - T_base) % T_incr
+ return sp.Piecewise(
+ (data_array[0], T < T_base),
+ (data_array[-1], T >= T_base + (len(data_array) - 1) * T_incr),
+ (((T_incr - mod) * data_symbol[idx] + (mod * data_symbol[idx + 1])) / T_incr, True)
+ )
+ elif parameter.func == 'lookup':
+ return density_lookup(T, *resolved_args)
+ elif parameter.func == 'thermal_diffusivity':
+ return thermal_diffusivity_func(*resolved_args)
+ else:
+ raise ValueError("No known format in the dict: parameter = ", parameter)
+ else:
+ return parameter # assume everything which is not a dict is a float, int, etc.
+
+
+# Material Data Class
+@dataclass
+class Ti6Al4V: # 90% Ti, 6% Al, 4% V
+ latent_heat = 290000 # m²/s²
+ temperature_solidus_1 = 1876.0 # k Ref [1]
+ temperature_solidus_2 = 1880.0 # k Ref [2]
+ temperature_liquidus = 1928.0 # k Ref [2]
+ density_room = 4430.0 # kg/m³ Ref [5]
+ thermal_expansion_coefficient = 10e-6 # /k Ref[6]
+ heat_capacity_solidus = 700.0 # J/kg-k
+ heat_capacity_liquidus = 1126.0 # J/kg-k
+ heat_conductivity_solidus = 21.58 # W/m-k
+ heat_conductivity_liquidus = 33.60 # W/m-k
+ density_temp_array = [1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0]
+ density_data_array = [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0]
+ density_temp_base = 1878.0
+ density_temp_incr = 8.0
+ # thermal_diffusivity_solidus = 727.7
+ # thermal_diffusivity_liquidus = 708.6
+ # electrical_resistivity = [1, 2, 5, 10, 20, 50, 100]
+ # electrical_resistivity_temperatures = [800, 900, 1200, 1500, 1700, 1900, 2000]
+ # Define temperature and density arrays
+ density_temp_array_lookup = np.array([1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0])
+ density_data_array_lookup = np.array([4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0])
+
+
+@dataclass
+class Constants:
+ temperature_room = 298.15 # k
+ N_a = 6.022141e23 # /mol
+ u = 1.660538e-27 # kg
+
+
+# Instantiate the Ti6Al4V class
+ti6al4v = Ti6Al4V()
+
+Ti6Al4V_latent_heat = ti6al4v.latent_heat
+Ti6Al4V_thermal_expansion_coefficient = ti6al4v.thermal_expansion_coefficient
+
+
+@dataclass
+class Solidus_Liquidus_Parameter:
+ sol: float
+ liq: float
+
+
+@dataclass
+class Functional_Parameter:
+ func: str
+ args: list
+
+
+# Parameters for Ti6Al4V
+Ti6Al4V_temperatures = Solidus_Liquidus_Parameter(
+ sol=0.5 * (ti6al4v.temperature_solidus_1 + ti6al4v.temperature_solidus_2),
+ liq=ti6al4v.temperature_liquidus)
+
+Ti6Al4V_heat_conductivity = Solidus_Liquidus_Parameter(
+ sol=ti6al4v.heat_conductivity_solidus,
+ liq=ti6al4v.heat_conductivity_liquidus)
+
+Ti6Al4V_specific_heat_capacity = Solidus_Liquidus_Parameter(
+ sol=ti6al4v.heat_capacity_solidus,
+ liq=ti6al4v.heat_capacity_liquidus)
+
+Ti6Al4V_density_tec = Functional_Parameter(
+ func='thermal_expansion',
+ args=[Constants.temperature_room, ti6al4v.density_room, ti6al4v.thermal_expansion_coefficient])
+
+Ti6Al4V_density_int = Functional_Parameter(
+ func='interpolation',
+ args=[ti6al4v.density_temp_array, ti6al4v.density_data_array])
+
+Ti6Al4V_density_in_array = Functional_Parameter(
+ func='interpolation_in_array', # TODO
+ args=[ti6al4v.density_temp_base, ti6al4v.density_temp_incr, ti6al4v.density_data_array,
+ sp.IndexedBase(ps.TypedSymbol('density_data_symbol', PointerType(sp.Basic('float32'))), # PointerType(create_type(np.float32)) PointerType(base_type=float32) PointerType(sp.Basic('float32'))
+ shape=len(ti6al4v.density_data_array))])
+
+Ti6Al4V_density_lookup = Functional_Parameter(
+ func='lookup',
+ args=[ti6al4v.density_temp_array_lookup, ti6al4v.density_data_array_lookup])
+
+Ti6Al4V_thermal_diffusivity = Functional_Parameter(
+ func='thermal_diffusivity',
+ args=[Ti6Al4V_heat_conductivity, Ti6Al4V_density_lookup, Ti6Al4V_specific_heat_capacity]
+ # solidus=ti6al4v.thermal_diffusivity_solidus,
+ # liquidus=ti6al4v.thermal_diffusivity_liquidus
+)
+'''
with CodeGeneration() as ctx:
data_type = "float64" if ctx.double_accuracy else "float32"
u, u_tmp = ps.fields(f"u, u_tmp: {data_type}[2D]", layout='fzyx')
- kappa = sp.Symbol("kappa")
+ thermal_diffusivity = sp.Symbol("thermal_diffusivity")
+ thermal_diffusivity_out = ps.fields(f"thermal_diffusivity_out: {data_type}[2D]", layout='fzyx')
dx = sp.Symbol("dx")
dt = sp.Symbol("dt")
- heat_pde = ps.fd.transient(u) - kappa * (ps.fd.diff(u, 0, 0) + ps.fd.diff(u, 1, 1))
+ heat_pde = ps.fd.transient(u) - thermal_diffusivity * (ps.fd.diff(u, 0, 0) + ps.fd.diff(u, 1, 1))
discretize = ps.fd.Discretization2ndOrder(dx=dx, dt=dt)
heat_pde_discretized = discretize(heat_pde)
heat_pde_discretized = heat_pde_discretized.args[1] + heat_pde_discretized.args[0].simplify()
+ thermal_diffusivity_expr = get_parameters(u.center(), Ti6Al4V_parameters.Ti6Al4V_temperatures, Ti6Al4V_parameters.Ti6Al4V_thermal_diffusivity)
- @ps.kernel
- def update():
- u_tmp.center @= heat_pde_discretized
-
-
- ac = ps.AssignmentCollection(update)
+ ac = ps.AssignmentCollection(
+ subexpressions=[
+ ps.Assignment(thermal_diffusivity, thermal_diffusivity_expr)
+ ],
+ main_assignments=[
+ ps.Assignment(u_tmp.center(), heat_pde_discretized),
+ ps.Assignment(thermal_diffusivity_out.center(), thermal_diffusivity)
+ ])
ac = ps.simp.simplifications.add_subexpressions_for_divisions(ac)
+ # ac = ps.simp.simplifications.add_subexpressions_for_field_reads(ac)
+ # ac = ps.simp.sympy_cse(ac)
- generate_sweep(ctx, 'HeatEquationKernel', ac)
+ generate_sweep(ctx, 'HeatEquationKernel', ac, varying_parameters=((data_type, str(thermal_diffusivity)),))