bugfix 2d example

Author: pancetta

pySDC/implementations/problem_classes/HeatEquation_2D_FD_periodic.py

@@ -2,7 +2,7 @@
 
 import numpy as np
 import scipy.sparse as sp
-from scipy.sparse.linalg import splu
+from scipy.sparse.linalg import cg
 
 from pySDC.core.Problem import ptype
 from pySDC.core.Errors import ParameterError, ProblemError
@@ -75,6 +75,9 @@ def __get_A(N, nu, dx):
         doffsets = np.concatenate((offsets, np.delete(offsets, zero_pos - 1) - N[0]))
 
         A = sp.diags(dstencil, doffsets, shape=(N[0], N[0]), format='csc')
+        # stencil = [1, -2, 1]
+        # A = sp.diags(stencil, [-1, 0, 1], shape=(N[0], N[0]), format='csc')
+
         A = sp.kron(A, sp.eye(N[0])) + sp.kron(sp.eye(N[1]), A)
         A *= nu / (dx ** 2)
 
@@ -93,7 +96,7 @@ def eval_f(self, u, t):
         """
 
         f = self.dtype_f(self.init)
-        f.values = self.A.dot(u.values)
+        f.values = self.A.dot(u.values.flatten()).reshape(self.params.nvars)
         return f
 
     def solve_system(self, rhs, factor, u0, t):
@@ -111,8 +114,9 @@ def solve_system(self, rhs, factor, u0, t):
         """
 
         me = self.dtype_u(self.init)
-        L = splu(sp.eye(self.params.nvars, format='csc') - factor * self.A)
-        me.values = L.solve(rhs.values)
+        me.values = cg(sp.eye(self.params.nvars[0] * self.params.nvars[1], format='csc') - factor * self.A,
+                       rhs.values.flatten(), x0=u0.values.flatten(), tol=1E-12)[0]
+        me.values = me.values.reshape(self.params.nvars)
         return me
 
     def u_exact(self, t):
@@ -128,7 +132,7 @@ def u_exact(self, t):
 
         me = self.dtype_u(self.init)
         xvalues = np.array([i * self.dx for i in range(self.params.nvars[0])])
-        me.values = np.kron(np.sin(np.pi * self.params.freq * xvalues), np.sin(np.pi * self.params.freq * xvalues)) * \
+        xv, yv = np.meshgrid(xvalues, xvalues)
+        me.values = np.sin(np.pi * self.params.freq * xv) * np.sin(np.pi * self.params.freq * yv) * \
             np.exp(-t * self.params.nu * (np.pi * self.params.freq) ** 2)
-        me.values = me.values.flatten()
         return me

pySDC/playgrounds/order_test/heat1d.py

@@ -0,0 +1,130 @@
+# import pySDC.helpers.plot_helper as plt_helper
+#
+# import pickle
+# import os
+import numpy as np
+
+from pySDC.implementations.datatype_classes.mesh import mesh, rhs_imex_mesh
+from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
+from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
+from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right
+from pySDC.implementations.controller_classes.allinclusive_multigrid_nonMPI import allinclusive_multigrid_nonMPI
+from pySDC.implementations.problem_classes.HeatEquation_1D_FD_periodic import heat1d_periodic
+from pySDC.implementations.problem_classes.HeatEquation_1D_FD_forced import heat1d_forced
+from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh
+
+from pySDC.helpers.stats_helper import filter_stats, sort_stats
+
+
+def main():
+    # initialize level parameters
+    level_params = dict()
+    level_params['restol'] = 1E-12
+    level_params['dt'] = None
+
+    # This comes as read-in for the step class (this is optional!)
+    step_params = dict()
+    step_params['maxiter'] = None
+
+    # This comes as read-in for the problem class
+    problem_params = dict()
+    problem_params['nu'] = 1.0
+    problem_params['freq'] = 2
+    problem_params['nvars'] = [2 ** 14 - 1]#, 2 ** 13]
+
+    # 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['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
+
+    # initialize controller parameters
+    controller_params = dict()
+    controller_params['logger_level'] = 30
+
+    # Fill description dictionary for easy hierarchy creation
+    description = dict()
+    description['problem_class'] = heat1d_forced
+    description['dtype_u'] = mesh
+    description['dtype_f'] = rhs_imex_mesh
+    description['sweeper_class'] = 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
+
+
+    # setup parameters "in time"
+    t0 = 0
+    Tend = 2.0
+
+    dt_list = [Tend / 2 ** i for i in range(0, 4)]
+    niter_list = [100]#[1, 2, 3, 4]
+
+    for niter in niter_list:
+
+        err = 0
+        for dt in dt_list:
+
+            print('Working with dt = %s and k = %s iterations...' % (dt, niter))
+
+            description['step_params']['maxiter'] = niter
+            description['level_params']['dt'] = dt
+
+            # instantiate the controller
+            controller = allinclusive_multigrid_nonMPI(num_procs=1, controller_params=controller_params,
+                                                       description=description)
+
+            # get initial values on finest level
+            P = controller.MS[0].levels[0].prob
+            uinit = P.u_exact(t0)
+
+            # call main function to get things done...
+            uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
+
+            # compute exact solution and compare
+            uex = P.u_exact(Tend)
+            err_new = abs(uex - uend)
+
+            print('   error at time %s: %s' % (Tend, err_new))
+            if err > 0:
+                print('   order of accuracy: %6.4f' % (np.log(err / err_new) / np.log(2)))
+
+            err = err_new
+
+            # # filter statistics by type (number of iterations)
+            # filtered_stats = filter_stats(stats, type='niter')
+            #
+            # # convert filtered statistics to list of iterations count, sorted by process
+            # iter_counts = sort_stats(filtered_stats, sortby='time')
+            #
+            # # compute and print statistics
+            # niters = np.array([item[1] for item in iter_counts])
+            # out = '   Mean number of iterations: %4.2f' % np.mean(niters)
+            # print(out)
+            # out = '   Range of values for number of iterations: %2i ' % np.ptp(niters)
+            # # f.write(out + '\n')
+            # print(out)
+            # out = '   Position of max/min number of iterations: %2i -- %2i' % \
+            #       (int(np.argmax(niters)), int(np.argmin(niters)))
+            # # f.write(out + '\n')
+            # print(out)
+            # out = '   Std and var for number of iterations: %4.2f -- %4.2f' % \
+            #       (float(np.std(niters)), float(np.var(niters)))
+            # # f.write(out + '\n')
+            # # f.write(out + '\n')
+            # print(out)
+        print()
+
+if __name__ == "__main__":
+    main()