From e6f8d15a614d42e0b7c6f6ef3590b5886702aa7d 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: Sun, 12 Mar 2023 00:13:51 +0100
Subject: [PATCH] Added 'solve_exactly()' for Burgers class.
---
Snakefile | 7 ++-
scripts/tcd/Equation.py | 102 +++++++++++++++++++++++++++++++++++++++-
2 files changed, 106 insertions(+), 3 deletions(-)
diff --git a/Snakefile b/Snakefile
index 6616d90..c4dceb2 100644
--- a/Snakefile
+++ b/Snakefile
@@ -23,18 +23,23 @@ TODO: Contemplate allowing vector input for ICs
TODO: Discuss how wavelet details should be plotted
TODO: Ask whether we are dealing with inviscid Burgers only
TODO: Ask whether it is helpful to enforce boundary at every SSPRK3 step
+ -> Done (should be not needed as we enforce for RHS)
TODO: Ask why we limit the initial time step
TODO: Ask whether cfl number should be absolute value
TODO: Ask why cells 1 and -2 are adjusted for Dirichlet boundary
+TODO: Ask about introducing subclasses of Burgers -> Done (do)
Urgent:
TODO: Add Burgers class completely
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 'solve_exactly()' to Burgers -> Done
TODO: Add 'update_right_hand_side()' to Burgers -> Done
TODO: Add Dirichlet boundary -> Done
+TODO: Clean up exact solution for Burgers
+TODO: Add subclasses of Burgers
+TODO: Ignore ghost domains by flagging all cells with it for extreme outlier
TODO: Add functions to create each object from dict
TODO: Initialize boundary config properly
TODO: Test Burgers
diff --git a/scripts/tcd/Equation.py b/scripts/tcd/Equation.py
index 5298b04..6b170e2 100644
--- a/scripts/tcd/Equation.py
+++ b/scripts/tcd/Equation.py
@@ -460,8 +460,8 @@ class Burgers(Equation):
"""
pass
- # num_periods = np.floor(self.wave_speed * self.final_time /
- # mesh.interval_len)
+ num_periods = np.floor(self.wave_speed * self.final_time /
+ mesh.interval_len)
#
# grid = np.repeat(mesh.non_ghost_cells, self._quadrature.num_nodes) + \
# mesh.cell_len / 2 * np.tile(self._quadrature.nodes, mesh.num_cells)
@@ -483,6 +483,104 @@ class Burgers(Equation):
# exact = np.reshape(exact, (1, exact.size))
#
# return grid, exact
+ grid = []
+ exact = []
+
+ for cell in range(len(mesh)):
+ eval_points = mesh[cell] + cell_len / 2 * self._quadrature.get_eval_points()
+ grid.append(eval_points)
+ eval_values = []
+ if scheme_label == 'Linear':
+ for point in range(len(eval_points)):
+ new_entry = self._init_cond.calculate(eval_points[point] - self._wave_speed * self._final_time
+ + num_periods * self._interval_len)
+ eval_values.append(new_entry)
+ exact.append(eval_values)
+
+ grid = np.reshape(np.array(grid), (1, len(grid) * len(grid[0])))
+
+ if scheme_label == 'Linear':
+ # u(x,t) = u(x - wavespeed*time, 0)
+ exact = np.reshape(np.array(exact), (1, len(exact) * len(exact[0])))
+
+ if scheme_label == 'Burgers' and initial_condition == 'Sine':
+ # u(x,t) = u(x - u0*time, 0)
+ exact = self._init_cond._height_adjustment + \
+ self._init_cond._stretch_factor * \
+ self.implicit_burgers_solver(
+ grid, self._final_time)
+ elif scheme_label == 'Burgers' and initial_condition == 'DiscontinuousConstant':
+ # u(x,t) = u(x - u0*time, 0)
+ exact = self.rarefaction_wave(grid, self._final_time)
+
+ def rarefaction_wave(self, grid_values, final_time):
+ uexact = 0 * grid_values
+ N = np.size(grid_values)
+
+ for i in range(N):
+ if grid_values[0, i] <= - final_time:
+ uexact[0, i] = -1
+ elif - final_time < grid_values[0, i] < final_time:
+ uexact[0, i] = grid_values[0, i] / final_time
+ else:
+ uexact[0, i] = 1
+
+ return uexact
+
+ def implicit_burgers_solver(self, grid_values, final_time):
+
+ ## Code Adapted from Hesthaven DG Code
+
+ sspeed = self._init_cond._height_adjustment
+ uexact = np.zeros(np.size(grid_values))
+
+ # initialize values
+ xx = np.linspace(0, 1, 200)
+ uu = np.sin(self._init_cond._factor * np.pi * xx)
+
+ te = final_time * (2 / self._interval_len)
+ xt0 = (2 / self._interval_len) * (
+ grid_values-sspeed * final_time).reshape(
+ np.size(grid_values))
+ for ix in range(len(xt0)):
+ if xt0[ix] > 1:
+ xt0[ix] -= 2 * np.floor(0.5 * (xt0[ix]+1))
+ elif xt0[ix] < -1:
+ xt0[ix] += 2 * np.floor(0.5 * (1-xt0[ix]))
+
+ ay = abs(xt0[ix])
+ # Finding the closest characteristic
+ i0 = 0
+ for j in range(200):
+ xt = xx[j]+uu[j] * te
+ if ay < xt:
+ break
+ else:
+ i0 = j
+
+ # This is our initial guess
+ us = uu[i0]
+ un = us
+ for k in range(400):
+ us = un
+ x0 = ay-us * te
+ un = us-(us-np.sin(self._init_cond._factor * np.pi * x0)) / (
+ 1+self._init_cond._factor * np.pi * np.cos(
+ self._init_cond._factor * np.pi * x0) * te)
+
+ if abs(un-us) < 1e-15:
+ break
+
+ y = np.sign(xt0[ix]) * un
+ if abs(ay-1) < 1e-15:
+ y = 0
+ if abs(final_time) < 1e-15:
+ y = np.sin(self._init_cond._factor * np.pi * xt0[ix])
+ uexact[ix] = y
+
+ burgers_exact = uexact.reshape((1, np.size(grid_values)))
+
+ return burgers_exact
@enforce_boundary()
def update_right_hand_side(self, projection: ndarray) -> ndarray:
--
GitLab