added problem definition and driver file for multiscale example

Author: Daniel Ruprecht

examples/acoustic_1d_imex/ProblemClass_multiscale.py

@@ -0,0 +1,174 @@
+r"""
+  One-dimensional IMEX acoustic-advection
+  =========================
+  
+  Integrate the linear 1D acoustic-advection problem:
+  
+  .. math::
+  u_t + U u_x + c p_x & = 0 \\
+  p_t + U p_x + c u_x & = 0.
+  
+"""
+
+import numpy as np
+import scipy.sparse as sp
+import scipy.sparse.linalg as LA
+
+from pySDC.Problem import ptype
+from pySDC.datatype_classes.mesh import mesh, rhs_imex_mesh
+
+from buildWave1DMatrix import getWave1DMatrix, getWave1DAdvectionMatrix
+
+class acoustic_1d_imex(ptype):
+    """
+    Example implementing the forced 1D heat equation with Dirichlet-0 BC in [0,1]
+
+    Attributes:
+      solver: Sharpclaw solver
+      state:  Sharclaw state
+      domain: Sharpclaw domain
+    """
+
+    def __init__(self, cparams, dtype_u, dtype_f):
+        """
+        Initialization routine
+
+        Args:
+            cparams: 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
+        assert 'nvars' in cparams
+        assert 'cs' in cparams
+        assert 'cadv' in cparams
+        assert 'order_adv' in cparams
+        assert 'multiscale' in cparams
+        
+        # add parameters as attributes for further reference
+        for k,v in cparams.items():
+            setattr(self,k,v)
+
+        # invoke super init, passing number of dofs, dtype_u and dtype_f
+        super(acoustic_1d_imex,self).__init__(self.nvars,dtype_u,dtype_f)
+        
+        self.mesh   = np.linspace(0.0, 1.0, self.nvars[1], endpoint=False)
+        self.dx     = self.mesh[1] - self.mesh[0]
+        
+        self.Dx     = -self.cadv*getWave1DAdvectionMatrix(self.nvars[1], self.dx, self.order_adv)
+        self.Id, A  = getWave1DMatrix(self.nvars[1], self.dx, ['periodic','periodic'], ['periodic','periodic'])
+        self.A      = -self.cs*A
+                
+    def solve_system(self,rhs,factor,u0,t):
+        """
+        Simple linear solver for (I-dtA)u = rhs
+
+        Args:
+            rhs: right-hand side for the nonlinear system
+            factor: abbrev. for the node-to-node stepsize (or any other factor required)
+            u0: initial guess for the iterative solver (not used here so far)
+            t: current time (e.g. for time-dependent BCs)
+
+        Returns:
+            solution as mesh
+        """
+        
+        M  = self.Id - factor*self.A
+        
+        b = np.concatenate( (rhs.values[0,:], rhs.values[1,:]) )
+        
+        sol = LA.spsolve(M, b)
+
+        me = mesh(self.nvars)
+        me.values[0,:], me.values[1,:] = np.split(sol, 2)
+        
+        return me
+
+
+    def __eval_fexpl(self,u,t):
+        """
+        Helper routine to evaluate the explicit part of the RHS
+
+        Args:
+            u: current values (not used here)
+            t: current time
+
+        Returns:
+            explicit part of RHS
+        """
+
+
+        b = np.concatenate( (u.values[0,:], u.values[1,:]) )
+        sol = self.Dx.dot(b)
+        
+        fexpl        = mesh(self.nvars)
+        fexpl.values[0,:], fexpl.values[1,:] = np.split(sol, 2)
+
+        return fexpl
+
+
+    def __eval_fimpl(self,u,t):
+        """
+        Helper routine to evaluate the implicit part of the RHS
+
+        Args:
+            u: current values
+            t: current time (not used here)
+
+        Returns:
+            implicit part of RHS
+        """
+
+        b = np.concatenate( (u.values[0,:], u.values[1,:]) )
+        sol = self.A.dot(b)
+        
+        fimpl             = mesh(self.nvars,val=0)
+        fimpl.values[0,:], fimpl.values[1,:] = np.split(sol, 2)
+        
+        return fimpl
+
+
+    def eval_f(self,u,t):
+        """
+        Routine to evaluate both parts of the RHS
+
+        Args:
+            u: current values
+            t: current time
+
+        Returns:
+            the RHS divided into two parts
+        """
+
+        f = rhs_imex_mesh(self.nvars)
+        f.impl = self.__eval_fimpl(u,t)
+        f.expl = self.__eval_fexpl(u,t)
+        return f
+
+    def u_exact(self,t):
+        """
+        Routine to compute the exact solution at time t
+
+        Args:
+            t: current time
+
+        Returns:
+            exact solution
+        """
+        
+        sigma_0 = 0.1
+        k       = 7.0*2.0*np.pi
+        x_0     = 0.75
+        x_1     = 0.25
+        
+        ms = 0.0
+        if self.multiscale:
+          ms = 1.0
+
+        me             = mesh(self.nvars)
+        me.values[0,:] = np.exp(-np.square(self.mesh-x_0-self.cs*t)/(sigma_0*sigma_0)) + ms*np.exp(-np.square(self.mesh-x_1-self.cs*t)/(sigma_0*sigma_0))*np.cos(k*(self.mesh-self.cs*t)/sigma_0)
+        me.values[1,:] = me.values[0,:]
+        return me
+
+

examples/acoustic_1d_imex/runmultiscale.py

@@ -0,0 +1,97 @@
+
+from pySDC import CollocationClasses as collclass
+
+import numpy as np
+
+from ProblemClass_multiscale import acoustic_1d_imex
+from examples.acoustic_1d_imex.HookClass import plot_solution
+
+from pySDC.datatype_classes.mesh import mesh, rhs_imex_mesh
+from pySDC.sweeper_classes.imex_1st_order import imex_1st_order
+import pySDC.PFASST_stepwise as mp
+from pySDC import Log
+from pySDC.Stats import grep_stats, sort_stats
+
+# Sharpclaw imports
+#from clawpack import pyclaw
+#from clawpack import riemann
+from matplotlib import pyplot as plt
+
+
+if __name__ == "__main__":
+
+    # set global logger (remove this if you do not want the output at all)
+    logger = Log.setup_custom_logger('root')
+
+    num_procs = 1
+
+    # This comes as read-in for the level class
+    lparams = {}
+    lparams['restol'] = 1E-10
+
+    sparams = {}
+    sparams['maxiter'] = 4
+
+    # setup parameters "in time"
+    t0   = 0.0
+    Tend = 3.0
+    dt = Tend/float(154)
+    
+    # This comes as read-in for the problem class
+    pparams = {}
+    pparams['nvars']      = [(2,512)]
+    pparams['cadv']       = 0.05
+    pparams['cs']         = 1.0
+    pparams['order_adv']  = 5
+    pparams['multiscale'] = True
+
+    # This comes as read-in for the transfer operations
+    tparams = {}
+    tparams['finter'] = True
+
+    # Fill description dictionary for easy hierarchy creation
+    description = {}
+    description['problem_class']     = acoustic_1d_imex
+    description['problem_params']    = pparams
+    description['dtype_u']           = mesh
+    description['dtype_f']           = rhs_imex_mesh
+    description['collocation_class'] = collclass.CollGaussLobatto
+    description['num_nodes']         = 3
+    description['sweeper_class']     = imex_1st_order
+    description['level_params']      = lparams
+    description['hook_class']        = plot_solution
+    #description['transfer_class'] = mesh_to_mesh_1d
+    #description['transfer_params'] = tparams
+
+    # quickly generate block of steps
+    MS = mp.generate_steps(num_procs,sparams,description)
+
+    # get initial values on finest level
+    P = MS[0].levels[0].prob
+    uinit = P.u_exact(t0)
+
+    # call main function to get things done...
+    uend,stats = mp.run_pfasst(MS,u0=uinit,t0=t0,dt=dt,Tend=Tend)
+
+    # compute exact solution and compare
+    uex = P.u_exact(Tend)
+
+    fig = plt.figure(figsize=(8,8))
+
+    sigma_0 = 0.1
+    k       = 7.0*2.0*np.pi
+    x_0     = 0.75
+    x_1     = 0.25
+
+    #plt.plot(P.mesh, uex.values[0,:],  '+', color='b', label='u (exact)')
+    plt.plot(P.mesh, uend.values[1,:], '-', color='b', label='SDC')
+    #plt.plot(P.mesh, uex.values[1,:],  '+', color='r', label='p (exact)')
+    #plt.plot(P.mesh, uend.values[1,:], '-', color='b', linewidth=2.0, label='p (SDC)')
+    p_slow = np.exp(-np.square(P.mesh-x_0-pparams['cadv']*Tend)/(sigma_0*sigma_0))
+    plt.plot(P.mesh, p_slow, '-', color='r', markersize=4, label='slow mode')
+    plt.legend(loc=2)
+    plt.xlim([0, 1])
+    plt.ylim([-0.1, 1.1])
+    fig.gca().grid()
+    plt.show()
+    #plt.gcf().savefig('fwsw-sdc-K'+str(sparams['maxiter'])+'-M'+str(description['num_nodes'])+'.pdf', bbox_inches='tight')