Coverage for pyrc \ validation \ coupling \ analytic.py: 0%
96 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-13 16:59 +0200
1# -------------------------------------------------------------------------------
2# Copyright (C) 2026 Joel Kimmich, Tim Jourdan
3# ------------------------------------------------------------------------------
4# License
5# This file is part of PyRC, distributed under GPL-3.0-or-later.
6# ------------------------------------------------------------------------------
8# Script for validation
9#
10# Calculation of the dynamic solution of two coupled volumes with different temperature
11#
12# Analytic solution
13#
14import numpy as np
16from pyrc.visualization.plot import LinePlot
17from pyrc.validation.coupling import setup
18from pyrc.validation.network import Block, Model
20thermal_conductivity = 10
21density = 10
22heat_capacity = 10
23a = thermal_conductivity / (density * heat_capacity)
26class ModelCoupling(Model):
27 def __init__(self, blocks: list[Block]):
28 """
29 Initialize heat equation model with blocks.
31 Parameters
32 ----------
33 blocks : list of Block
34 List of Block objects defining initial conditions
35 """
36 super().__init__(blocks)
38 def get_function(self, n_terms=60):
39 """
40 Returns function u(x,t) for 1D heat equation with adiabatic boundaries.
42 Parameters
43 ----------
44 n_terms : int
45 Number of Fourier series terms
47 Returns
48 -------
49 callable
50 Function u(x, t) that returns temperature for given positions x and times t.
51 The function returns a value, vector or matrix.
52 Same x values are in rows, same t values are in columns. (x, t)
53 """
54 L = self.length
56 # Calculate A0
57 A0 = sum(self.widths * np.array([b.initial_temperature for b in self.blocks])) / L
59 # Calculate An coefficients
60 def calculate_An(n):
61 A_n = 0
62 k_n = n * np.pi / L
63 for j, block in enumerate(self.blocks):
64 T_j = block.initial_temperature
65 x1 = self.x_positions[j]
66 x2 = self.x_positions[j + 1]
67 A_n += T_j * (np.sin(k_n * x2) - np.sin(k_n * x1))
69 return (2 / (L * k_n)) * A_n
71 # Precompute coefficients
72 coeffs = [calculate_An(n) for n in range(1, n_terms + 2)]
74 def u(x, t):
75 x = np.atleast_1d(x)
76 t = np.atleast_1d(t)
77 X, T = np.meshgrid(x, t, indexing="ij")
79 result = np.full_like(X, A0)
81 for n in range(1, n_terms + 2):
82 result += coeffs[n - 1] * np.cos(n * np.pi * X / L) * np.exp(-self.alpha * (n * np.pi / L) ** 2 * T)
84 # Override t=0 with exact initial conditions
85 t_zero_indices = np.where(t == 0)[0]
86 if len(t_zero_indices) > 0:
87 for t_idx in t_zero_indices:
88 for i, block in enumerate(self.blocks):
89 if i == len(self.blocks) - 1: # Last block includes right endpoint
90 mask = (x >= self.x_positions[i]) & (x <= self.x_positions[i + 1])
91 else:
92 mask = (x >= self.x_positions[i]) & (x < self.x_positions[i + 1])
93 result[mask, t_idx] = block.initial_temperature
95 return result.squeeze()
97 return u
99 def get_function_new(self, n_terms=100):
100 """
101 Returns function u(x,t) for 1D heat equation with adiabatic boundaries.
103 Parameters
104 ----------
105 n_terms : int
106 Number of Fourier series terms
108 Returns
109 -------
110 callable
111 Function u(x, t) that returns temperature for given positions x and times t.
112 The function returns a value, vector or matrix.
113 Same x values are in rows, same t values are in columns. (x, t)
114 """
116 L = self.length
118 def sum_constant(k):
119 """
120 Calculates constant values for each block for a given number of increment.
122 Parameters
123 ----------
124 k : int
125 The increment number.
127 Returns
128 -------
129 list :
130 A list with the constants for all blocks for the given icrement number k.
131 """
132 result = []
134 for n, block in enumerate(self.blocks):
135 T_n = block.initial_temperature
136 x_n = self.x_positions[n]
137 x_n_plus_one = self.x_positions[n + 1]
139 result.append(
140 T_n * 2 / (np.pi * k) * (np.sin(k * np.pi * x_n_plus_one / L) - np.sin(k * np.pi * x_n / L))
141 )
142 return result
144 def temperature_constants():
145 # Counter of blocks with constant initial temperautres start at 0 here, not from 1 like in the derivation
146 # from Juri Joussen
147 result = []
148 for n, block in enumerate(self.blocks):
149 result.append(block.width / L * block.initial_temperature)
150 return result
152 temperature_const = temperature_constants()
153 sum_constants = [sum_constant(n) for n in range(1, n_terms + 1)]
155 def u(x, t):
156 x = np.atleast_1d(x)
157 t = np.atleast_1d(t)
158 X, T = np.meshgrid(x, t, indexing="ij")
160 result = np.full_like(X, np.array(temperature_const))
162 for k in range(1, n_terms + 2):
163 # The material properties are missing?!
164 result += sum_constants[k - 1] * np.cos(k * np.pi * X / L) * np.exp(-((k * np.pi / L) ** 2) * T)
166 # Override t=0 with exact initial conditions
167 t_zero_indices = np.where(t == 0)[0]
168 if len(t_zero_indices) > 0:
169 for t_idx in t_zero_indices:
170 for i, block in enumerate(self.blocks):
171 if i == len(self.blocks) - 1: # Last block includes right endpoint
172 mask = (x >= self.x_positions[i]) & (x <= self.x_positions[i + 1])
173 else:
174 mask = (x >= self.x_positions[i]) & (x < self.x_positions[i + 1])
175 result[mask, t_idx] = block.initial_temperature
177 return result.squeeze()
179 return u
181 def error(self, n_terms, time):
182 return np.exp(-self.alpha * (n_terms * np.pi / 15) ** 2 * time)
184 def get_n_terms(self, time, maximum_error=1e-10):
185 return int((-np.log(maximum_error) / (self.alpha * time)) ** 0.5 * self.length / np.pi + 1)
188def run_coupling_analytic(plot_data: str = None, n_terms=None):
189 blocks, _, times, name_add_on = setup()
191 model = ModelCoupling(blocks)
193 if n_terms is None:
194 n_terms = model.get_n_terms(time=0.1, maximum_error=1e-12)
195 else:
196 n_terms = int(n_terms)
197 function = model.get_function(n_terms=n_terms)
199 result = function(x=model.x_values, t=np.array(times))
201 if plot_data:
202 plot = LinePlot(
203 x=model.x_values,
204 ys=result.T,
205 labels=[f"{float(t):.1f} s" for t in times],
206 y_title="Temperature / °C (analytic)",
207 x_title="Length / m",
208 )
209 plot.plot()
210 plot.ax.set_xlim(left=0, right=model.x_positions[-1])
211 plot.fig.legend(loc="outside upper right", ncols=5)
212 plot.save(plot_data)
213 return model.x_values, result.T, times, name_add_on
216if __name__ == "__main__":
217 run_cooupling_analytic(True)