From 520bb6f25fe09be8cbe0f6aa9281f757a8ad2916 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?K=C3=BChle=2C=20Laura=20Christine=20=28lakue103=29?=
<laura.kuehle@uni-duesseldorf.de>
Date: Fri, 10 Mar 2023 15:43:36 +0100
Subject: [PATCH] Added 'update_right_hand_side()' for Burgers class.
---
Snakefile | 6 +-
scripts/tcd/Equation.py | 138 ++++++++++++++++++++++++++++++++++++++++
2 files changed, 142 insertions(+), 2 deletions(-)
diff --git a/Snakefile b/Snakefile
index b415bad..03a1175 100644
--- a/Snakefile
+++ b/Snakefile
@@ -32,9 +32,11 @@ TODO: Add 'update_time_step()' to Burgers -> Done
TODO: Add derivative basis to Basis class -> Done
TODO: Fix wrong type for number of quadrature nodes -> Done
TODO: Add 'solve_exactly()' to Burgers
-TODO: Add 'update_right_hand_side()' to Burgers
-TODO: Add initialization to Burgers
+TODO: Add 'update_right_hand_side()' to Burgers -> Done
TODO: Add Dirichlet boundary
+TODO: Test Burgers
+TODO: Rework Burgers for speed-up and better structure
+TODO: Add initialization to Burgers
TODO: Use cfl_number for updating, not just time (equation-related?)
Critical, but not urgent:
diff --git a/scripts/tcd/Equation.py b/scripts/tcd/Equation.py
index d0bd5ce..5298b04 100644
--- a/scripts/tcd/Equation.py
+++ b/scripts/tcd/Equation.py
@@ -7,6 +7,7 @@ from typing import Tuple
from abc import ABC, abstractmethod
import numpy as np
from numpy import ndarray
+from sympy import Symbol
from .Basis_Function import Basis
from .Boundary_Condition import enforce_boundary
@@ -16,6 +17,9 @@ from .projection_utils import do_initial_projection
from .Quadrature import Quadrature
+x = Symbol('x')
+
+
class Equation(ABC):
"""Abstract class for equation.
@@ -496,6 +500,23 @@ class Burgers(Equation):
"""
pass
+ volume_integral_matrix = self._burgers_volume_integral(projection)
+ boundary_matrix = self._burgers_boundary_matrix(projection)
+
+ # Initialize vector and set first entry to accommodate for ghost cell
+ right_hand_side = [0]
+
+ for j in range(self._mesh.num_cells):
+ right_hand_side.append(2*(volume_integral_matrix[:, j + 1] +
+ boundary_matrix[:, j + 1]))
+
+ # Set ghost cells to respective value
+ # (Periodic, Updated to Dirichlet in enforce_boundary_condition
+ # for DiscontinuousConstant problem)
+ # right_hand_side[0] = right_hand_side[self._mesh.num_cells]
+ # right_hand_side.append(right_hand_side[1])
+
+ return np.transpose(right_hand_side)
# right_hand_side = np.zeros_like(projection)
# if self._wave_speed > 0:
# right_hand_side[:, self._mesh.num_ghost_cells:
@@ -516,3 +537,120 @@ class Burgers(Equation):
# -self._mesh.num_ghost_cells])
#
# return right_hand_side
+
+ def _burgers_volume_integral(self, projection):
+ basis_vector = self._basis.basis
+ derivative_basis_vector = self._basis.derivative
+
+ # Initialize matrix and set first entry to accommodate for ghost cell
+ volume_integral_matrix = [0]
+
+ for cell in range(self._mesh.num_cells):
+ new_row = []
+ # Approximation u at Gauss-Legendre points
+ approx_eval = [sum(projection[degree][cell + 1]
+ * basis_vector[degree].subs(
+ x, self._quadrature.nodes[point])
+ for degree in range(
+ self._basis.polynomial_degree + 1))
+ for point in range(self._quadrature.nodes)]
+ approx_eval = np.array(approx_eval)
+
+ # print("u", approx_eval)
+ # Quadrature Evaluation
+ for degree in range(self._basis.polynomial_degree + 1):
+ weighted_derivatives = [
+ derivative_basis_vector[degree].subs(
+ x, self._quadrature.nodes[point])
+ * self._quadrature.weights[point]
+ for point in range(self._quadrature.num_nodes)]
+ new_entry = sum(list(approx_eval**2 * weighted_derivatives)) / 2
+
+ new_row.append(np.float64(new_entry))
+
+ new_row = np.array(new_row).T
+ volume_integral_matrix.append(new_row)
+
+ # Set ghost cells to respective value
+ # (Periodic, Updated to Dirichlet in enforce_boundary_condition
+ # for DiscontinuousConstant problem)
+ # volume_integral_matrix[0] = volume_integral_matrix[
+ # self._mesh.num_cells]
+ # volume_integral_matrix.append(volume_integral_matrix[1])
+
+ return np.transpose(np.array(volume_integral_matrix))
+
+ def _burgers_local_Lax_Friedrichs(self, projection):
+ # Calculating the left and right boundaries for each cell,
+ # u_{j-1/2}^+, u_{j+1/2}^-
+ boundary_left, boundary_right = self._calculate_boundary_points(
+ projection, self._basis.polynomial_degree + 1)
+ # print("shape BL", boundary_left.shape)
+
+ # Initializing Burgers local Lax-Friedrichs flux matrix
+ burgers_flux = [0]
+ for j in range(self._mesh.num_cells):
+ # approximations j+1/2^-, j+1/2^+ for interior cells
+ # and max velocity
+ approx_minus = boundary_right[:, j+1]
+ approx_plus = boundary_left[:, j + 2]
+ max_velocity = max(abs(approx_minus), abs(approx_plus))
+
+ # respective fluxes
+ flux_minus = approx_minus**2 / 2
+ flux_plus = approx_plus**2 / 2
+
+ # local Lax-Friedrichs flux
+ numerical_flux = 0.5 * (flux_minus + flux_plus - max_velocity *
+ (approx_plus - approx_minus))
+
+ burgers_flux.append(numerical_flux)
+
+ # Set Ghost cells to respective value
+ # (Periodic, Updated to Dirichlet in enforce_boundary_condition
+ # for DiscontinuousConstant problem)
+ # burgers_flux[0] = burgers_flux[self._mesh.num_cells]
+ # burgers_flux.append(burgers_flux[1])
+
+ return np.array(burgers_flux)
+
+ def _burgers_boundary_matrix(self, projection):
+ burgers_LF_flux = self._burgers_local_Lax_Friedrichs(projection)
+
+ # Initializing boundary matrix
+ boundary_matrix = [0]
+
+ for j in range(self._mesh.num_cells):
+ new_row = []
+ for degree in range(self._basis.polynomial_degree + 1):
+ new_row.append(np.sqrt(degree + 0.5) *
+ (-burgers_LF_flux[j + 1] + (-1)**degree *
+ burgers_LF_flux[j]))
+
+ boundary_matrix.append(new_row)
+
+ # Set Ghost cells to respective value
+ # (Periodic, Updated to Dirichlet in enforce_boundary_condition
+ # for DiscontinuousConstant problem)
+ # boundary_matrix[0] = boundary_matrix[self._mesh.num_cells]
+ # boundary_matrix.append(boundary_matrix[1])
+
+ boundary_matrix = np.transpose(np.array(boundary_matrix))
+
+ return boundary_matrix[0]
+
+ @staticmethod
+ def _calculate_boundary_points(projection, max_degree):
+ # Approximation at j-1/2 ^+
+ boundary_left = [sum(projection[degree][cell] * (-1) ** degree *
+ np.sqrt(degree + 0.5)
+ for degree in range(max_degree))
+ for cell in range(len(projection[0]))]
+
+ # Approximation at j+1/2 ^-
+ boundary_right = [sum(projection[degree][cell]* np.sqrt(degree + 0.5)
+ for degree in range(max_degree))
+ for cell in range(len(projection[0]))]
+
+ return np.reshape(np.array(boundary_left), (1, len(boundary_left))), \
+ np.reshape(np.array(boundary_right), (1, len(boundary_right)))
--
GitLab