using print(...) instead of print ..., for python3, added some comments

Author: BenjaminRodenberg

coupling_schemes.py

@@ -15,7 +15,13 @@
 
 def estimate_coupling_neumann_BC(left_domain, u_left, right_domain, u_right):
     """
-    @type left_domain Domain
+    estimate neumann boundary condition at the coupling interface between left_domain and right_domain from existing nodal data
+    :param left_domain: left domain
+    :type left_domain: Domain
+    :param u_left: solution on left domain
+    :param left_domain: right domain
+    :type right_domain: Domain
+    :param u_right: solution on right domain
     @type right_domain Domain
     """
     # set neumann BC at coupling interface
@@ -27,7 +33,7 @@ def estimate_coupling_neumann_BC(left_domain, u_left, right_domain, u_right):
             u_neumann_coupled__ = finite_difference_schemes.central_FD_at(left_domain.u.shape[0]-1, right_domain.grid.h, u_glued, order=numeric_parameters.neumann_coupling_order)
         else: # use modified operator for non-identical mesh size
             if numeric_parameters.neumann_coupling_order != 2:
-                print "Operator of order %d is not implemented!!!" % numeric_parameters.neumann_coupling_order
+                print("Operator of order %d is not implemented!!!" % numeric_parameters.neumann_coupling_order)
                 quit()
             fraction = left_domain.grid.h / right_domain.grid.h  # normalize to right domain's meshwidth
             p = np.array([-fraction, 0, 1.0])
@@ -36,13 +42,17 @@ def estimate_coupling_neumann_BC(left_domain, u_left, right_domain, u_right):
             u = np.array([u_left[-2], u_right[0], u_right[1]])
             u_neumann_coupled__ = 1.0/right_domain.grid.h * (u.dot(c))
     else:
-        print "not implemented schemes for coupling Neumann BC demanded!"
+        print("not implemented schemes for coupling Neumann BC demanded!")
         quit()
 
     return u_neumann_coupled__
 
 
 class CouplingScheme(object):
+    """
+    abstract class defining a coupling scheme.
+    """
+
     __metaclass__ = abc.ABCMeta
 
     name = "Coupling Scheme"
@@ -53,23 +63,28 @@ def __init__(self):
 
     @abc.abstractmethod
     def perform(self, t0, tau, left_domain, right_domain):
+        """
+        abstract method to perform one coupled timestep with coupling between equation on left_domain and right_domain
+        :param t0: current time
+        :param tau: time step size
+        :param left_domain: left domain with independent model
+        :param right_domain: right domain with independent model
+        :return:
+        """
         return
 
 
 class FullyImplicitCoupling(CouplingScheme):
+    """
+    fully implicit coupling with iterative coupling
+    """
 
     name = "Fully Implicit Coupling"
 
     def __init__(self):
         super(FullyImplicitCoupling, self).__init__()
 
     def perform(self, t0, tau, left_domain, right_domain):
-        """
-        @type left_domain Domain
-        @type right_domain Domain
-        :param t0:
-        """
-
         residual = np.inf
         tol = numeric_parameters.fixed_point_tol
         n_steps_max = numeric_parameters.n_max_fixed_point_iterations
@@ -120,7 +135,7 @@ def perform(self, t0, tau, left_domain, right_domain):
             n_steps += 1
 
         if n_steps == n_steps_max:
-            print "maximum number of steps exceeded!"
+            print("maximum number of steps exceeded!")
             return False
 
         # update solution
@@ -133,6 +148,9 @@ def perform(self, t0, tau, left_domain, right_domain):
 
 
 class FullyExplicitCoupling(CouplingScheme):
+    """
+    fully explicit coupling with staggered approach
+    """
 
     name = "Fully Explicit Coupling"
 
@@ -195,6 +213,9 @@ def perform(self, t0, tau, left_domain, right_domain):
 
 
 class WaveformCoupling(CouplingScheme):
+    """
+    waveform coupling relying on waveform relaxation
+    """
 
     name = "Waveform Coupling"
 
@@ -295,7 +316,7 @@ def perform(self, t0, tau, left_domain, right_domain):
             u_right_new[0] = u_dirichlet_continuous_m1(t0+tau)  # we have to set the (known and changing) dirichlet value manually, since this value is not changed by the timestepping
 
             if numeric_parameters.neumann_coupling_order != 2:
-                print "Operator of order %d is not implemented!!!" % numeric_parameters.neumann_coupling_order
+                print("Operator of order %d is not implemented!!!" % numeric_parameters.neumann_coupling_order)
                 quit()
 
             fraction = left_domain.grid.h / right_domain.grid.h  # normalize to right domain's meshwidth
@@ -320,6 +341,9 @@ def perform(self, t0, tau, left_domain, right_domain):
 
 
 class ExplicitPredictorCoupling(CouplingScheme):
+    """
+    predictor coupling using a predictor scheme for the coupled participants
+    """
 
     name = "Explicit Predictor Coupling"
 
@@ -393,6 +417,9 @@ def perform(self, t0, tau, left_domain, right_domain):
 
 
 class SemiImplicitExplicitCoupling(CouplingScheme):
+    """
+    coupling scheme using a combination of explicit and implicit coupling
+    """
 
     name = "Semi Implicit Explicit Coupling"
 
@@ -470,7 +497,7 @@ def perform(self, t0, tau, left_domain, right_domain):
             n_steps += 1
 
         if n_steps == n_steps_max:
-            print "maximum number of steps exceeded!"
+            print("maximum number of steps exceeded!")
             return False
 
         # update solution
@@ -483,6 +510,9 @@ def perform(self, t0, tau, left_domain, right_domain):
 
 
 class StrangSplittingCoupling(CouplingScheme):
+    """
+    coupling using Strang splitting
+    """
 
     name = "Strang Splitting Coupling"
 
@@ -565,7 +595,9 @@ def perform(self, t0, tau, left_domain, right_domain):
 
 
 class MonolithicScheme(CouplingScheme):
-    # todo this class hierarchy stinks!
+    """
+    monolithic scheme for solving an equation without coupling
+    """
 
     name = "Monolithic Approach"
 

finite_difference_schemes.py

@@ -15,9 +15,11 @@ def central_FD_at(i, h, u, derivative=1, order=2):
         if order == 2:
             du = 1.0/(2.0 * h) * (-1.0 * u[i-1] + 1.0 * u[i+1])
         else:
-            print "not implemented!"
+            print("central_FD_at(...) with derivarive == 1 and order != 2 not implemented!")
+            quit()
     else:
-        print "not implemented!"
+        print("central_FD_at(...) with derivarive > 1 not implemented!")
+        quit()
 
     return du
 
@@ -51,8 +53,10 @@ def one_sided_forward_FD_at(i, h, u, derivative=1, order=1):
         elif order == 5:
             du = 1.0/h * (-137.0/60.0 * u[i] + 5.0 * u[i+1] - 5.0 * u[i+2] + 10.0/3.0 * u[i+3] - 5.0/4.0 * u[i+4] + 1.0/5.0 * u[i+5])
         else:
-            print "not implemented!"
+            print("one_sided_forward_FD_at(...) with derivarive == 1 and order > 5 not implemented!")
+            quit()
     else:
-        print "not implemented!"
+        print("one_sided_forward_FD_at(...) with derivarive > 1 not implemented!")
+        quit()
 
     return du

main.py

@@ -38,36 +38,36 @@ def experimental_series(taus, u_ref_left, u_ref_right, experiment):
     """
     @type taus list
     """
-    print "EXPERIMENTAL SERIES"
+    print("EXPERIMENTAL SERIES")
 
     errors_left = []
     errors_right = []
     experiment_counter = 1
     n_experiments = taus.__len__()
 
     for tau in taus:
-        print "----"
-        print "%i of %i" % (experiment_counter, n_experiments)
-        print "tau = %f" % tau
+        print("----")
+        print("%i of %i" % (experiment_counter, n_experiments))
+        print("tau = %f" % tau)
 
         u_left, u_right = experiment(tau)
 
         # compute error
         e_tot_left = np.sum(np.abs(u_left-u_ref_left)) / u_left.shape[0]
         e_tot_right = np.sum(np.abs(u_right-u_ref_right)) / (u_right.shape[0])
 
-        print "Error total = %.2e | %.2e" % (e_tot_left, e_tot_right)
+        print("Error total = %.2e | %.2e" % (e_tot_left, e_tot_right))
 
         experiment_counter += 1
         errors_left += [e_tot_left]
         errors_right += [e_tot_right]
 
-    print "----"
+    print("----")
     np.set_printoptions(formatter={'float': lambda x: format(x, '6.3E')}) # print in scientific notation
-    print repr(np.array(errors_left))
-    print repr(np.array(errors_right))
-    print taus
-    print
+    print(repr(np.array(errors_left)))
+    print(repr(np.array(errors_right)))
+    print(taus)
+    print()
 
     return errors_left, errors_right
 
@@ -145,7 +145,7 @@ def solve_monolithic_problem(tau, time_stepping_scheme, domain):
 def compute_reference_solution(domain, tau, time_integration_scheme = time_integration.ImplicitTrapezoidalRule()):
     # compute monolithic at constant spatial and very fine temporal resolution
 
-    print "REFERENCE SOLUTION for grid "+str(domain.grid.x)
+    print("REFERENCE SOLUTION for grid "+str(domain.grid.x))
 
     return solve_monolithic_problem(tau, time_integration_scheme, domain)
 
@@ -166,7 +166,7 @@ def coupling_experiment(coupling_scheme, time_stepping_schemes, left_domain, rig
     """
     if type(time_stepping_schemes) is list:
         experiment_name = time_stepping_schemes[0].name + " - " + time_stepping_schemes[1].name + " - " + coupling_scheme.name
-        print experiment_name
+        print(experiment_name)
         experiment_errors_left, experimental_errors_right = experimental_series(experiment_timesteps, u_ref_left, u_ref_right, lambda tau: solve_inhomogeneously_coupled_problem(
             tau, time_stepping_schemes, coupling_scheme, [left_domain, right_domain]))
     else:
@@ -175,7 +175,7 @@ def coupling_experiment(coupling_scheme, time_stepping_schemes, left_domain, rig
             experiment_name += " - "+coupling_scheme.name_suffix
         except AttributeError:
             pass
-        print experiment_name
+        print(experiment_name)
         experiment_errors_left, experimental_errors_right = experimental_series(experiment_timesteps, u_ref_left, u_ref_right, lambda tau: solve_coupled_problem(tau, time_stepping_schemes, coupling_scheme, left_domain, right_domain))
 
     return Experiment(experiment_name, [left_domain.grid.h, right_domain.grid.h], experiment_timesteps, experiment_errors_left, experimental_errors_right, additional_numerical_parameters=numeric_parameters.numeric_parameters_dict)
@@ -193,7 +193,7 @@ def monolithic_experiment(time_stepping_scheme, monolithic_domain, u_ref_left, u
     :return:
     """
     experiment_name = time_stepping_scheme.name + " - " + coupling_schemes.MonolithicScheme().name
-    print experiment_name
+    print(experiment_name)
     experiment_errors_left, experimental_errors_right = experimental_series(experiment_timesteps, u_ref_left, u_ref_right, lambda tau: solve_monolithic_problem(tau, time_stepping_scheme, monolithic_domain))
 
     if type(monolithic_domain.grid) is IrregularGrid:
@@ -250,7 +250,7 @@ def save_experiments(experimental_series, prefix):
 u_ref_left_regular, u_ref_right_regular = compute_reference_solution(monolithic_domain_regular, tau_ref, time_integration_scheme=time_integration.RungeKutta4)
 u_ref_left_nonregular, u_ref_right_nonregular = compute_reference_solution(monolithic_domain_nonregular, tau_ref, time_integration_scheme=time_integration.RungeKutta4)
 
-print "### experimental series 1: order degradation with classical schemes on regular domain"
+print("### experimental series 1: order degradation with classical schemes on regular domain")
 
 experiments = list()
 experiments.append(monolithic_experiment(time_integration.ExplicitHeun, monolithic_domain_regular, u_ref_left_regular, u_ref_right_regular, experiment_timesteps))
@@ -261,7 +261,7 @@ def save_experiments(experimental_series, prefix):
 experiments.append(coupling_experiment(coupling_schemes.FullyImplicitCoupling(), time_integration.RungeKutta4, left_domain, right_domain_coarse, u_ref_left_regular, u_ref_right_regular, experiment_timesteps))
 save_experiments(experiments, 'series1_Order')
 
-print "### experimental series 2: Customized schemes Semi Implicit-Explicit coupling and Predictor coupling on regular domain"
+print("### experimental series 2: Customized schemes Semi Implicit-Explicit coupling and Predictor coupling on regular domain")
 
 experiments = list()
 experiments.append(monolithic_experiment(time_integration.ExplicitHeun, monolithic_domain_regular, u_ref_left_regular, u_ref_right_regular, experiment_timesteps))
@@ -272,7 +272,7 @@ def save_experiments(experimental_series, prefix):
 experiments.append(coupling_experiment(coupling_schemes.SemiImplicitExplicitCoupling(), time_integration.ImplicitTrapezoidalRule, left_domain, right_domain_coarse, u_ref_left_regular, u_ref_right_regular, experiment_timesteps))
 save_experiments(experiments, 'series2_Custom')
 
-print "### experimental series 3: Strang splitting coupling on non-regular domain"
+print("### experimental series 3: Strang splitting coupling on non-regular domain")
 
 experiments = list()
 experiments.append(monolithic_experiment(time_integration.ExplicitHeun, monolithic_domain_nonregular, u_ref_left_nonregular, u_ref_right_nonregular, experiment_timesteps))
@@ -284,7 +284,7 @@ def save_experiments(experimental_series, prefix):
 experiments.append(coupling_experiment(coupling_schemes.StrangSplittingCoupling(), [time_integration.RungeKutta4, time_integration.ImplicitTrapezoidalRule], left_domain, right_domain_fine, u_ref_left_nonregular, u_ref_right_nonregular, experiment_timesteps))
 save_experiments(experiments, 'series3_Strang')
 
-print "### experimental series 4: Waveform relaxation coupling on non-regular domain"
+print("### experimental series 4: Waveform relaxation coupling on non-regular domain")
 
 experiments = list()
 experiments.append(monolithic_experiment(time_integration.ImplicitTrapezoidalRule, monolithic_domain_nonregular, u_ref_left_nonregular, u_ref_right_nonregular, experiment_timesteps))

postproc.py

@@ -17,7 +17,8 @@
 name = parser.parse_args().foldername
 plottingdir = r'./' + name
 filenames = glob.glob(plottingdir+"/*.json")
-print filenames
+print("reading from")
+print(filenames)
 
 marker_collection = ('o','x','d','s','^','<','>','v','.')
 
@@ -27,9 +28,6 @@
 if plot_right:
     marker_right = itertools.cycle(marker_collection)
 
-for filename in filenames:
-    print filename
-
 if plot_right:
     f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12,9))
 else:
@@ -58,8 +56,8 @@ def concat_w_zero_padding(a,b):
 
 
 for filename in filenames:
-    print "---"
-    print filename
+    print("---")
+    print(filename)
     with open(filename) as data_file:
 
         data = json.load(data_file)

space_discretization.py

@@ -43,7 +43,7 @@ def second_derivative_matrix(grid, **kwargs):
     elif type(grid) is IrregularGrid:
         return irregular_matrix_operator(grid.x, **kwargs)
     else:
-        print "unsupported grid type!"
+        print("in second_derivative_matrix(...) unsupported grid type!")
         quit()
 
 
@@ -83,7 +83,8 @@ def irregular_matrix_operator(x, **kwargs):
         # manually add it to R
         R[1] = 1.0/h**2 * kwargs.get('dirichlet_l')
     else:
-        print "insufficient boundary conditions at left boundary!"
+        print("in irregular_matrix_operator(...): insufficient boundary conditions at left boundary!")
+        quit()
 
     if 'neumann_r' in kwargs:
         # modify row and column corresponding to neumann BC at u_N
@@ -99,7 +100,8 @@ def irregular_matrix_operator(x, **kwargs):
         # manually add it to R
         R[-2] = 1.0/h**2 * kwargs.get('dirichlet_r')
     else:
-        print "insufficient boundary conditions at right boundary!"
+        print("in irregular_matrix_operator(...): insufficient boundary conditions at right boundary!")
+        quit()
 
     return A, R
 
@@ -132,7 +134,8 @@ def regular_matrix_operator(h, N, **kwargs):
         # manually add it to R
         R[1] = 1.0/h**2 * kwargs.get('dirichlet_l')
     else:
-        print "insufficient boundary conditions at left boundary!"
+        print("in regular_matrix_operator(...): insufficient boundary conditions at left boundary!")
+        quit()
 
     if 'neumann_r' in kwargs:
         # modify row and column corresponding to neumann BC at u_N
@@ -146,6 +149,7 @@ def regular_matrix_operator(h, N, **kwargs):
         # manually add it to R
         R[-2] = 1.0/h**2 * kwargs.get('dirichlet_r')
     else:
-        print "insufficient boundary conditions at right boundary!"
+        print("in regular_matrix_operator(...): insufficient boundary conditions at right boundary!")
+        quit()
 
     return A, R