added FPUT

Author: pancetta

pySDC/implementations/problem_classes/FermiPastaUlamTsingou.py

@@ -0,0 +1,122 @@
+from __future__ import division
+from numba import jit
+import numpy as np
+
+from pySDC.core.Problem import ptype
+from pySDC.implementations.datatype_classes.particles import particles, acceleration
+from pySDC.core.Errors import ParameterError
+
+
+# noinspection PyUnusedLocal
+class fermi_pasta_ulam_tsingou(ptype):
+    """
+    Example implementing the outer solar system problem
+    """
+
+    def __init__(self, problem_params, dtype_u, dtype_f):
+        """
+        Initialization routine
+
+        Args:
+            problem_params (dict): custom parameters for the example
+            dtype_u: particle data type (will be passed parent class)
+            dtype_f: acceleration data type (will be passed parent class)
+        """
+
+        # these parameters will be used later, so assert their existence
+        essential_keys = ['npart', 'alpha', 'k', 'energy_modes']
+        for key in essential_keys:
+            if key not in problem_params:
+                msg = 'need %s to instantiate problem, only got %s' % (key, str(problem_params.keys()))
+                raise ParameterError(msg)
+
+        # invoke super init, passing nparts, dtype_u and dtype_f
+        super(fermi_pasta_ulam_tsingou, self).__init__(problem_params['npart'], dtype_u, dtype_f, problem_params)
+
+    @staticmethod
+    @jit
+    def fast_acceleration(N, alpha, pos, accel):
+        for n in range(1, N - 1):
+            accel[n - 1] = pos[n - 1] - 2.0 * pos[n] + pos[n + 1] + \
+                alpha * ((pos[n + 1] - pos[n]) ** 2 - (pos[n] - pos[n - 1]) ** 2)
+
+    def eval_f(self, u, t):
+        """
+        Routine to compute the RHS
+
+        Args:
+            u (dtype_u): the particles
+            t (float): current time (not used here)
+        Returns:
+            dtype_f: RHS
+        """
+        me = acceleration(self.init, val=0.0)
+
+        self.fast_acceleration(self.params.npart, self.params.alpha, u.pos.values, me.values[1:self.params.npart - 1])
+
+        return me
+
+    def u_exact(self, t):
+        """
+        Routine to compute the exact/initial trajectory at time t
+
+        Args:
+            t (float): current time
+        Returns:
+            dtype_u: exact/initial position and velocity
+        """
+        assert t == 0.0, 'error, u_exact only works for the initial time t0=0'
+
+        me = particles(self.init)
+
+        for n in range(self.params.npart):
+            me.pos.values[n] = np.sin(self.params.k * np.pi * n / (self.params.npart - 1))
+            me.vel.values[n] = 0.0
+
+        return me
+
+    @staticmethod
+    @jit
+    def fast_hamiltonian(N, alpha, pos, vel):
+        ham = 0.0
+        for n in range(N - 1):
+            ham += 0.5 * vel[n] ** 2 + 0.5 * (pos[n + 1] - pos[n]) ** 2 + alpha / 3.0 * (pos[n + 1] - pos[n]) ** 3
+        return ham
+
+    def eval_hamiltonian(self, u):
+        """
+        Routine to compute the Hamiltonian
+
+        Args:
+            u (dtype_u): the particles
+        Returns:
+            float: hamiltonian
+        """
+        return self.fast_hamiltonian(self.params.npart, self.params.alpha, u.pos.values, u.vel.values)
+
+    def eval_mode_energy(self, u):
+        """
+        Routine to compute the energy following
+        http://www.scholarpedia.org/article/Fermi-Pasta-Ulam_nonlinear_lattice_oscillations
+
+        Args:
+            u (dtype_u): the particles
+        Returns:
+            dict: energies
+        """
+
+        energy = {}
+
+        for k in self.params.energy_modes:
+
+            Qk = np.sqrt(2.0 / (self.params.npart - 1)) * \
+                sum([u.pos.values[n] * np.sin(np.pi * k * n / (self.params.npart - 1))
+                     for n in range(1, self.params.npart)])
+            Qkdot = np.sqrt(2.0 / (self.params.npart - 1)) * \
+                sum([u.vel.values[n] * np.sin(np.pi * k * n / (self.params.npart - 1))
+                     for n in range(1, self.params.npart)])
+
+            omegak2 = 4.0 * np.sin(k * np.pi / (2.0 * (self.params.npart - 1))) ** 2
+            energy[k] = 0.5 * (Qkdot ** 2 + omegak2 * Qk ** 2)
+
+        return energy

pySDC/projects/Hamiltonian/fput.py

@@ -0,0 +1,250 @@
+from __future__ import division
+from collections import defaultdict
+import os
+import dill
+import numpy as np
+
+import pySDC.helpers.plot_helper as plt_helper
+
+from pySDC.implementations.collocation_classes.gauss_lobatto import CollGaussLobatto
+from pySDC.implementations.controller_classes.allinclusive_classic_nonMPI import allinclusive_classic_nonMPI
+from pySDC.implementations.sweeper_classes.verlet import verlet
+from pySDC.implementations.datatype_classes.particles import particles, acceleration
+from pySDC.implementations.problem_classes.FermiPastaUlamTsingou import fermi_pasta_ulam_tsingou
+from pySDC.implementations.transfer_classes.TransferParticles_NoCoarse import particles_to_particles
+
+from pySDC.helpers.stats_helper import filter_stats, sort_stats
+from pySDC.projects.Hamiltonian.hamiltonian_and_energy_output import hamiltonian_and_energy_output
+
+
+def setup_fput():
+    """
+    Helper routine for setting up everything for the Fermi-Pasta-Ulam-Tsingou problem
+
+    Returns:
+        description (dict): description of the controller
+        controller_params (dict): controller parameters
+    """
+
+    # initialize level parameters
+    level_params = dict()
+    level_params['restol'] = 1E-08
+    level_params['dt'] = 2.5
+
+    # initialize sweeper parameters
+    sweeper_params = dict()
+    sweeper_params['collocation_class'] = CollGaussLobatto
+    sweeper_params['num_nodes'] = [5]
+    sweeper_params['spread'] = False
+
+    # initialize problem parameters for the Penning trap
+    problem_params = dict()
+    problem_params['npart'] = 34
+    problem_params['alpha'] = 0.25
+    problem_params['k'] = 1.0
+    problem_params['energy_modes'] = [[1, 2, 3, 4]]
+
+    # initialize step parameters
+    step_params = dict()
+    step_params['maxiter'] = 50
+
+    # initialize controller parameters
+    controller_params = dict()
+    controller_params['hook_class'] = hamiltonian_and_energy_output
+    controller_params['logger_level'] = 30
+    controller_params['predict'] = False
+
+    # Fill description dictionary for easy hierarchy creation
+    description = dict()
+    description['problem_class'] = fermi_pasta_ulam_tsingou
+    description['problem_params'] = problem_params
+    description['dtype_u'] = particles
+    description['dtype_f'] = acceleration
+    description['sweeper_class'] = verlet
+    description['sweeper_params'] = sweeper_params
+    description['level_params'] = level_params
+    description['step_params'] = step_params
+    description['space_transfer_class'] = particles_to_particles
+
+    return description, controller_params
+
+
+def run_simulation():
+    """
+    Routine to run the simulation of a second order problem
+
+    """
+
+    description, controller_params = setup_fput()
+    # set time parameters
+    t0 = 0.0
+    Tend = 5000.0
+    num_procs = 1
+
+    f = open('fput_out.txt', 'w')
+    out = 'Running fput problem with %s processors...' % num_procs
+    f.write(out + '\n')
+    print(out)
+
+    # instantiate the controller
+    controller = allinclusive_classic_nonMPI(num_procs=num_procs, controller_params=controller_params,
+                                             description=description)
+
+    # get initial values on finest level
+    P = controller.MS[0].levels[0].prob
+    uinit = P.u_exact(t=t0)
+
+    # call main function to get things done...
+    uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
+
+    # 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
+    for item in iter_counts:
+        out = 'Number of iterations for time %4.2f: %2i' % item
+        f.write(out)
+        print(out)
+
+    niters = np.array([item[1] for item in iter_counts])
+    out = '   Mean number of iterations: %4.2f' % np.mean(niters)
+    f.write(out + '\n')
+    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')
+    print(out)
+
+    # get runtime
+    timing_run = sort_stats(filter_stats(stats, type='timing_run'), sortby='time')[0][1]
+    out = '... took %s seconds to run this.' % timing_run
+    f.write(out + '\n')
+    print(out)
+    f.close()
+
+    # assert np.mean(niters) <= maxmeaniter, 'Mean number of iterations is too high, got %s' % np.mean(niters)
+
+    fname = 'data/fput.dat'
+    f = open(fname, 'wb')
+    dill.dump(stats, f)
+    f.close()
+
+    assert os.path.isfile(fname), 'Run for %s did not create stats file' % prob
+
+
+def show_results(cwd=''):
+    """
+    Helper function to plot the error of the Hamiltonian
+
+    Args:
+        cwd (str): current working directory
+    """
+
+    # read in the dill data
+    f = open(cwd + 'data/fput.dat', 'rb')
+    stats = dill.load(f)
+    f.close()
+
+    plt_helper.mpl.style.use('classic')
+    plt_helper.setup_mpl()
+
+    # HAMILTONIAN PLOTTING #
+    # extract error in hamiltonian and prepare for plotting
+    extract_stats = filter_stats(stats, type='err_hamiltonian')
+    result = defaultdict(list)
+    for k, v in extract_stats.items():
+        result[k.iter].append((k.time, v))
+    for k, v in result.items():
+        result[k] = sorted(result[k], key=lambda x: x[0])
+
+    plt_helper.newfig(textwidth=238.96, scale=0.89)
+
+    # Rearrange data for easy plotting
+    err_ham = 1
+    for k, v in result.items():
+        time = [item[0] for item in v]
+        ham = [item[1] for item in v]
+        err_ham = ham[-1]
+        plt_helper.plt.semilogy(time, ham, '-', lw=1, label='Iter ' + str(k))
+    print(err_ham)
+    # assert err_ham < 7.5E-15, 'Error in the Hamiltonian is too large for %s, got %s' % (prob, err_ham)
+
+    plt_helper.plt.xlabel('Time')
+    plt_helper.plt.ylabel('Error in Hamiltonian')
+    plt_helper.plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
+
+    fname = 'data/fput_hamiltonian'
+    plt_helper.savefig(fname)
+
+    assert os.path.isfile(fname + '.pdf'), 'ERROR: plotting did not create PDF file'
+    assert os.path.isfile(fname + '.pgf'), 'ERROR: plotting did not create PGF file'
+    assert os.path.isfile(fname + '.png'), 'ERROR: plotting did not create PNG file'
+
+    # ENERGY PLOTTING #
+    # extract error in hamiltonian and prepare for plotting
+    extract_stats = filter_stats(stats, type='energy_step')
+    result = sort_stats(extract_stats, sortby='time')
+
+    plt_helper.newfig(textwidth=238.96, scale=0.89)
+
+    # Rearrange data for easy plotting
+
+    for mode in result[0][1].keys():
+        time = [item[0] for item in result]
+        energy = [item[1][mode] for item in result]
+        plt_helper.plt.plot(time, energy, label=str(mode) + 'th mode')
+
+    plt_helper.plt.xlabel('Time')
+    plt_helper.plt.ylabel('Energy')
+    plt_helper.plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
+
+    fname = 'data/fput_energy'
+    plt_helper.savefig(fname)
+
+    assert os.path.isfile(fname + '.pdf'), 'ERROR: plotting did not create PDF file'
+    assert os.path.isfile(fname + '.pgf'), 'ERROR: plotting did not create PGF file'
+    assert os.path.isfile(fname + '.png'), 'ERROR: plotting did not create PNG file'
+
+    # POSITION PLOTTING #
+    # extract positions and prepare for plotting
+    extract_stats = filter_stats(stats, type='position')
+    result = sort_stats(extract_stats, sortby='time')
+
+    plt_helper.newfig(textwidth=238.96, scale=0.89)
+
+    # Rearrange data for easy plotting
+    nparts = len(result[0][1])
+    nsteps = len(result)
+    pos = np.zeros((nparts, nsteps))
+    time = np.zeros(nsteps)
+    for idx, item in enumerate(result):
+        time[idx] = item[0]
+        for n in range(nparts):
+            pos[n, idx] = item[1][n]
+
+    for n in range(nparts):
+        plt_helper.plt.plot(time, pos[n, :])
+
+    fname = 'data/fput_positions'
+    plt_helper.savefig(fname)
+
+    assert os.path.isfile(fname + '.pdf'), 'ERROR: plotting did not create PDF file'
+    assert os.path.isfile(fname + '.pgf'), 'ERROR: plotting did not create PGF file'
+    assert os.path.isfile(fname + '.png'), 'ERROR: plotting did not create PNG file'
+
+
+def main():
+    run_simulation()
+    show_results()
+
+if __name__ == "__main__":
+    main()