Added 'solve_exactly()' for Burgers class.

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

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: