Bug fix

aux_funcs.py

@@ -0,0 +1,40 @@
+def peclet_number(a, cell_widths, d):
+    """
+    Calculate the Peclet number for a given system.
+
+    Parameters
+    ----------
+    a : float
+        Advection velocity.
+    cell_widths : array_like
+        Array of cell widths.
+    d : float
+        Diffusion coefficient.
+
+    Returns
+    -------
+    ndarray
+        Array of Peclet numbers for each cell.
+    """
+    return a * cell_widths / d
+
+
+def CFL_condition(a, k, cell_widths):
+    """
+    Calculate the CFL condition for a given system.
+
+    Parameters
+    ----------
+    a : float
+        Advection velocity.
+    k : float
+        Time step factor.
+    cell_widths : array_like
+        Array of cell widths.
+
+    Returns
+    -------
+    ndarray
+        Array of CFL conditions for each cell.
+    """
+    return a * k / cell_widths

bcs.py

@@ -0,0 +1,76 @@
+from functools import partial
+
+# ROBIN BOUNDARY CONDITONS
+
+
+def robin_boundary_condition_matrix_elements_left(m, a, d):
+    # Left side Robin boundary coefficients for matrix equation
+    b1 = 1.0 / m.h(0) * (-a.p(0) * m.h(1) / (2 * m.hp(0)) - d.p(0) / m.hp(0))
+    c1 = 1.0 / m.h(0) * (-a.p(0) * m.h(0) / (2 * m.hp(0)) + d.p(0) / m.hp(0))
+
+    # Index and element value
+    locations = [(0, 0), (0, 1)]
+    values = (b1, c1)
+    return tuple([list(x) for x in zip(locations, values)])
+
+
+def robin_boundary_condition_matrix_elements_right(m, a, d):
+    # Right hand side Robin boundary coefficients for matrix equation
+    J = m.n_cells
+    aJ = (
+        1.0
+        / m.h(J - 1)
+        * (a.m(J - 1) * m.h(J - 1) / (2 * m.hm(J - 1)) + d.m(J - 1) / m.hm(J - 1))
+    )
+    bJ = (
+        1.0
+        / m.h(J - 1)
+        * (a.m(J - 1) * m.h(J - 2) / (2 * m.hm(J - 1)) - d.m(J - 1) / m.hm(J - 1))
+    )
+
+    J = m.n_cells  # Index and element value
+
+    # Index and element value
+    locations = [(J - 1, J - 2), (J - 1, J - 1)]
+    values = (aJ, bJ)
+    return tuple([list(x) for x in zip(locations, values)])
+
+
+def robin_boundary_condition_vector_elements(location, m, flux):
+    location = [location]
+    value = [flux / m.h(location)]
+    return tuple([list(x) for x in zip(location, value)])
+
+
+robin_boundary_condition_vector_elements_left = partial(
+    robin_boundary_condition_vector_elements, location=0
+)
+
+robin_boundary_condition_vector_elements_right = partial(
+    robin_boundary_condition_vector_elements, location=-1
+)
+
+# DIRICHLET BOUNDARY CONDITIONS
+
+
+def dirichlet_boundary_conditions(locations, values):
+    # Dirichlet boundary coefficients
+    return tuple([list(x) for x in zip(locations, values)])
+
+
+# Left hand side Dirichlet boundary coefficients for matrix equation
+dirichlet_boundary_condition_matrix_elements_left = partial(
+    dirichlet_boundary_conditions, locations=[(0, 0), (0, 1)], values=(0, 1)
+)
+
+dirichlet_boundary_condition_matrix_elements_right = partial(
+    dirichlet_boundary_conditions, locations=[(-1, -2), (-1, -1)], values=(0, 1)
+)
+
+dirichlet_boundary_condition_vector_elements_left = partial(
+    dirichlet_boundary_conditions, locations=(0,), values=(0,)
+)
+
+dirichlet_boundary_condition_vector_elements_right = partial(
+    dirichlet_boundary_conditions, locations=(-1,), values=(0,)
+)

cell_variable.py

@@ -0,0 +1,39 @@
+import numpy as np
+
+from .mesh import Mesh
+
+
+class CellVariable(np.ndarray):
+    """Representation of a variable defined at the cell centers. Provides interpolation functions to calculate the value at cell faces."""
+
+    def __new__(cls, input_array, mesh: Mesh):
+        # If `input_array` is actually just a constant
+        # convert it to an array of len the number of cells.
+        try:
+            len(input_array)
+        except TypeError:
+            input_array = input_array * np.ones(len(mesh.cells))
+
+        obj = np.asarray(input_array).view(cls)
+        obj.mesh = mesh
+        return obj
+
+    def __array_finalize__(self, obj: Mesh):
+        if obj is None:
+            return
+        self.mesh = getattr(obj, "mesh", None)
+        self.__get_items__ = getattr(obj, "__get_items__", None)
+
+    def m(self, i):
+        """Linear interpolation of the cell value at the right hand face i.e. along the _m_inus direction."""
+        return (
+            self.mesh.h(i) / (2 * self.mesh.hm(i)) * self[i - 1]
+            + self.mesh.h(i - 1) / (2 * self.mesh.hm(i)) * self[i]
+        )
+
+    def p(self, i):
+        """Linear interpolation of the cell value at the right hand face i.e. along the _p_lus direction."""
+        return (
+            self.mesh.h(i + 1) / (2 * self.mesh.hp(i)) * self[i]
+            + self.mesh.h(i) / (2 * self.mesh.hp(i)) * self[i + 1]
+        )

constant_flux.py

@@ -0,0 +1,94 @@
+import matplotlib
+
+# matplotlib.use("Agg")
+import matplotlib.animation as manimation
+import matplotlib.pyplot as plt
+import numpy as np
+from tqdm import tqdm
+
+from .cell_variable import CellVariable
+from .mesh import Mesh
+from .model import AdvectionDiffusionModel
+
+if __name__ == "__main__":
+    # https://www.mathworks.com/help/symbolic/heating-of-finite-slab.html
+    video = True
+
+    # L = 1.0  # m
+    # k = 10.0  # W/(m·K)
+    # rho = 1000.0  # kg/m³
+    # c = 1000.0  # J/(kg·K) # heat capacity
+    alpha = 117e-6  # k / (rho * c)  # Thermal diffusivity
+
+    mesh = Mesh.simple_grading(0, 0.6, n_points=50, ratio=0.1)
+    print(mesh.faces)
+
+    a = CellVariable(0, mesh=mesh)  # Advection velocity
+    d = CellVariable(alpha, mesh=mesh)  # Diffusion coefficient m2/s
+
+    time_step = 0.1e-1  # Time step, s
+    total_time = 120  # seconds
+    theta = 0  # 1 -> implicit
+    left_flux = +3e5 * alpha / 400  # W/m2
+    right_flux = 0
+
+    # Initial conditions
+    w_init = np.zeros_like(mesh.cells) + 20  # degrees
+
+    model = AdvectionDiffusionModel(
+        a,
+        d,
+        time_step,
+        mesh,
+        discretization="upwind",
+        theta=theta,
+    )
+    model.set_boundary_conditions(left_flux=left_flux, right_flux=right_flux)
+
+    model.initialize_system(w_init)
+
+    if not video:
+        w = model.solve_to_time(total_time)
+        fig, ax = plt.subplots()
+        ax.plot(mesh.cells, w_init, label="IC")
+        ax.plot(mesh.cells, w, label=f"t={total_time}", marker="o")
+        plt.show()
+    else:
+        FFMpegWriter = manimation.writers["ffmpeg"]
+        metadata = dict(
+            title="Movie Test", artist="Matplotlib", comment="Movie support!"
+        )
+        writer = FFMpegWriter(fps=15, metadata=metadata)
+
+        fig, ax = plt.subplots()
+        l1 = ax.plot([], [], "k-o", markersize=4)[0]
+        ann = ax.annotate(f"time 0", xy=(1, 1))
+
+        plt.xlim(np.min(mesh.faces), np.max(mesh.faces))
+        plt.ylim(0, 120)
+        l1.set_data(mesh.cells, w_init)
+
+        n_steps = int(total_time / time_step)
+        with writer.saving(fig, "fixed_heatflux.mp4", 100):
+            for i in tqdm(range(n_steps)):
+                w = model.take_step()
+
+                if i == 0:
+                    l1.set_data(mesh.cells, w_init)
+                    writer.grab_frame()
+
+                if i % 100 == 0 or i == 0:
+                    l1.set_data(mesh.cells, w)
+                    # l0.set_data(analytical_x, analytical_solution2)
+                    area = np.sum(w * mesh.cell_widths)
+
+                    if ann is not None:
+                        ann.remove()
+                    ann = ax.annotate(f"time {i*time_step}", xy=(1, 1))
+
+                    # print("#%d; t=%g; area=%g:" % (i, i * k, area))
+                    writer.grab_frame()
+
+        total_time = n_steps * time_step
+        print(w)
+        print(f"Simulated time: {total_time} s")

constant_temperature.py

@@ -0,0 +1,85 @@
+import matplotlib
+
+# matplotlib.use("Agg")
+import matplotlib.animation as manimation
+import matplotlib.pyplot as plt
+import numpy as np
+from tqdm import tqdm
+
+from .cell_variable import CellVariable
+from .mesh import Mesh
+from .model import AdvectionDiffusionModel
+
+if __name__ == "__main__":
+    # https://www.mathworks.com/help/symbolic/heating-of-finite-slab.html
+    video = False
+
+    mesh = Mesh.uniform(-1, 1, n_points=10)
+    print(mesh.faces)
+
+    a = CellVariable(0, mesh=mesh)  # Advection velocity
+    d = CellVariable(1, mesh=mesh)  # Diffusion coefficient m2/s
+
+    time_step = 1e-3  # Time step, s
+    total_time = 1  # seconds
+    theta = 1
+    left_value = 1
+    right_value = 1
+
+    # Initial conditions
+    w_init = np.zeros_like(mesh.cells)  # degrees
+    w_init[0] = left_value  # degrees
+    w_init[-1] = right_value  # degrees
+
+    model = AdvectionDiffusionModel(
+        a, d, time_step, mesh, discretization="upwind", theta=theta
+    )
+    model.set_boundary_conditions(left_value=left_value, right_value=right_value)
+    # model.set_boundary_conditions(left_flux=left_flux, right_flux=left_flux)
+
+    model.initialize_system(w_init)
+
+    if not video:
+        w = model.solve_to_time(total_time)
+        fig, ax = plt.subplots()
+        ax.plot(mesh.cells, w_init, label="IC")
+        ax.plot(mesh.cells, w, label=f"t={total_time}")
+        plt.show()
+
+    else:
+        FFMpegWriter = manimation.writers["ffmpeg"]
+        metadata = dict(title="Movie Test", artist="Matplotlib", comment="Movie support!")
+        writer = FFMpegWriter(fps=15, metadata=metadata)
+
+        fig, ax = plt.subplots()
+        l1 = ax.plot([], [], "k-o", markersize=4)[0]
+        ann = ax.annotate(f"time 0", xy=(1, 1))
+
+        plt.xlim(np.min(mesh.faces), np.max(mesh.faces))
+        plt.ylim(0, 1)
+        l1.set_data(mesh.cells, w_init)
+
+        n_steps = int(total_time / time_step)
+        with writer.saving(fig, "fixed_temperatures.mp4", 100):
+            for i in tqdm(range(n_steps)):
+                w = model.take_step()
+
+                if i == 0:
+                    l1.set_data(mesh.cells, w_init)
+                    writer.grab_frame()
+
+                if i % 1 == 0 or i == 0:
+                    l1.set_data(mesh.cells, w)
+                    # l0.set_data(analytical_x, analytical_solution2)
+                    area = np.sum(w * mesh.cell_widths)
+
+                    if ann is not None:
+                        ann.remove()
+                    ann = ax.annotate(f"time {i*time_step}", xy=(1, 1))
+
+                    # print("#%d; t=%g; area=%g:" % (i, i * k, area))
+                    writer.grab_frame()
+
+        total_time = n_steps * time_step
+        print(w)
+        print(f"Simulated time: {total_time} s")