bugfix...

Author: pancetta

pySDC/implementations/problem_classes/HeatEquation_1D_FD_periodic.py

@@ -125,8 +125,9 @@ def u_exact(self, t):
 
         if self.params.freq >= 0:
             xvalues = np.array([i * self.dx for i in range(self.params.nvars)])
+            rho = (2.0 - 2.0 * np.cos(np.pi * self.params.freq * self.dx)) / self.dx ** 2
             me.values = np.sin(np.pi * self.params.freq * xvalues) * \
-                np.exp(-t * self.params.nu * (np.pi * self.params.freq) ** 2)
+                np.exp(-t * self.params.nu * rho)
         else:
             np.random.seed(1)
             me.values = np.random.rand(self.params.nvars)

pySDC/playgrounds/order_test/heat1d.py

@@ -28,48 +28,49 @@ def main():
 
     # This comes as read-in for the problem class
     problem_params = dict()
-    problem_params['nu'] = 1.0
+    problem_params['nu'] = 0.01
     problem_params['freq'] = 2
-    problem_params['nvars'] = [2 ** 14 - 1]#, 2 ** 13]
+    problem_params['nvars'] = [2 ** 7]#, 2 ** 6]
 
     # This comes as read-in for the sweeper class
     sweeper_params = dict()
     sweeper_params['collocation_class'] = CollGaussRadau_Right
     sweeper_params['num_nodes'] = 3
     sweeper_params['QI'] = 'IE'
-    sweeper_params['spread'] = False
+    sweeper_params['spread'] = True
     sweeper_params['do_coll_update'] = False
 
     # initialize space transfer parameters
     space_transfer_params = dict()
     space_transfer_params['rorder'] = 2
     space_transfer_params['iorder'] = 2
-    space_transfer_params['periodic'] = False
+    space_transfer_params['periodic'] = True
 
     # initialize controller parameters
     controller_params = dict()
-    controller_params['logger_level'] = 30
+    controller_params['logger_level'] = 20
+    controller_params['predict'] = False
 
     # Fill description dictionary for easy hierarchy creation
     description = dict()
-    description['problem_class'] = heat1d_forced
+    description['problem_class'] = heat1d_periodic
     description['dtype_u'] = mesh
-    description['dtype_f'] = rhs_imex_mesh
-    description['sweeper_class'] = imex_1st_order
+    description['dtype_f'] = mesh#rhs_imex_mesh
+    description['sweeper_class'] = generic_implicit#imex_1st_order
     description['sweeper_params'] = sweeper_params
     description['step_params'] = step_params
     description['level_params'] = level_params
     description['problem_params'] = problem_params
-    # description['space_transfer_class'] = mesh_to_mesh  # pass spatial transfer class
-    # description['space_transfer_params'] = space_transfer_params  # pass paramters for spatial transfer
+    description['space_transfer_class'] = mesh_to_mesh  # pass spatial transfer class
+    description['space_transfer_params'] = space_transfer_params  # pass paramters for spatial transfer
 
 
     # setup parameters "in time"
-    t0 = 0
-    Tend = 2.0
+    t0 = 0.0
+    Tend = 1.0
 
-    dt_list = [Tend / 2 ** i for i in range(0, 4)]
-    niter_list = [100]#[1, 2, 3, 4]
+    dt_list = [Tend / 2 ** i for i in range(3, 6)]
+    niter_list = [2]#[1, 2, 3, 4]
 
     for niter in niter_list:
 
@@ -81,6 +82,8 @@ def main():
             description['step_params']['maxiter'] = niter
             description['level_params']['dt'] = dt
 
+            Tend = t0 + dt
+
             # instantiate the controller
             controller = allinclusive_multigrid_nonMPI(num_procs=1, controller_params=controller_params,
                                                        description=description)
@@ -93,6 +96,8 @@ def main():
             uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
 
             # compute exact solution and compare
+            # uex = compute_collocation_solution(controller)
+            # print(abs(uex - P.u_exact(Tend)))
             uex = P.u_exact(Tend)
             err_new = abs(uex - uend)
 
@@ -126,5 +131,25 @@ def main():
             # print(out)
         print()
 
+
+def compute_collocation_solution(controller):
+
+    Q = controller.MS[0].levels[0].sweep.coll.Qmat[1:, 1:]
+    A = controller.MS[0].levels[0].prob.A.todense()
+    dt = controller.MS[0].levels[0].dt
+
+    N = controller.MS[0].levels[0].prob.init
+    M = controller.MS[0].levels[0].sweep.coll.num_nodes
+
+    u0 = np.kron(np.ones(M), controller.MS[0].levels[0].u[0].values)
+
+    C = np.eye(M * N) - dt * np.kron(Q, A)
+
+    tmp = np.linalg.solve(C, u0)
+    print(np.linalg.norm(C.dot(tmp) - u0, np.inf))
+    uex = controller.MS[0].levels[0].prob.dtype_u(N)
+    uex.values[:] = tmp[(M-1) * N:]
+    return uex
+
 if __name__ == "__main__":
     main()

pySDC/tests/test_2d_fd_accuracy.py

@@ -59,10 +59,9 @@ def run_accuracy_check(nvars_list,problem_params):
         # create a mesh instance and fill it with the Laplacian of the sine wave
         u_lap = prob.dtype_u(init=prob.init)
         u_lap.values = -2*(np.pi * prob.params.freq) ** 2 * prob.params.nu * np.kron(np.sin(np.pi * prob.params.freq * xvalues), np.sin(np.pi * prob.params.freq * xvalues))
-        u_lap.values = u_lap.values.flatten()
 
         # compare analytic and computed solution using the eval_f routine of the problem class
-        err = abs(prob.eval_f(u,0) - u_lap)
+        err = abs(prob.eval_f(u, 0) - u_lap)
 
         # get id for this nvars and put error into dictionary
         id = ID(nvars=nvars)