Merge pull request #75 from amanabt/advection_1d_upwind_flux

Upwind flux implementation for 1D advection.

dg_maxwell/isoparam.py

@@ -42,6 +42,7 @@ def isoparam_x_2D(x_nodes, xi, eta):
     '''
     Finds the :math:`x` coordinate using isoparametric mapping of a
     :math:`2^{nd}` order element with :math:`8` nodes
+
     .. math:: (P_0, P_1, P_2, P_3, P_4, P_5, P_6, P_7)
 
     Here :math:`P_i` corresponds to :math:`(\\xi_i, \\eta_i)` coordinates,

dg_maxwell/lagrange.py

@@ -13,12 +13,12 @@
 
 def LGL_points(N):
     '''
-    Calculates : math: `N` Legendre-Gauss-Lobatto (LGL) points.
+    Calculates :math:`N` Legendre-Gauss-Lobatto (LGL) points.
     LGL points are the roots of the polynomial 
 
-    :math: `(1 - \\xi ** 2) P_{n - 1}'(\\xi) = 0`
+    .. math:: (1 - \\xi^2) P_{n - 1}'(\\xi) = 0
 
-    Where :math: `P_{n}(\\xi)` are the Legendre polynomials.
+    Where :math:`P_{n}(\\xi)` are the Legendre polynomials.
     This function finds the roots of the above polynomial.
 
     Parameters
@@ -85,11 +85,11 @@ def lobatto_weights(n):
 
 def gauss_nodes(n):
     '''
-    Calculates :math: `N` Gaussian nodes used for Integration by
+    Calculates :math:`N` Gaussian nodes used for Integration by
     Gaussia quadrature.
-    Gaussian node :math: `x_i` is the `i^{th}` root of
-    :math: `P_n(\\xi)`
-    Where :math: `P_{n}(\\xi)` are the Legendre polynomials.
+    Gaussian node :math:`x_i` is the :math:`i^{th}` root of
+    :math:`P_n(\\xi)`
+    Where :math:`P_{n}(\\xi)` are the Legendre polynomials.
 
     Parameters
     ----------
@@ -101,7 +101,7 @@ def gauss_nodes(n):
     -------
 
     gauss_nodes : numpy.ndarray
-                  The Gauss nodes :math: `x_i`.
+                  The Gauss nodes :math:`x_i`.
 
     **See:** A Wikipedia article about the Gauss-Legendre quadrature `here`_
     
@@ -163,21 +163,21 @@ def lagrange_polynomials(x):
     ----------
     
     x : numpy.array [N_LGL 1 1 1]
-        Contains the :math: `x` nodes using which the
+        Contains the :math:`x` nodes using which the
         lagrange basis functions need to be evaluated.
 
     Returns
     -------
     
     lagrange_basis_poly   : list
-                            A list of size `x.shape[0]` containing the
+                            A list of size ``x.shape[0]`` containing the
                             analytical form of the Lagrange basis polynomials
                             in numpy.poly1d form. This list is used in
                             integrate() function which requires the analytical
                             form of the integrand.
 
     lagrange_basis_coeffs : numpy.ndarray
-                            A :math: `N \\times N` matrix containing the
+                            A :math:`N \\times N` matrix containing the
                             coefficients of the Lagrange basis polynomials such
                             that :math:`i^{th}` lagrange polynomial will be the
                             :math:`i^{th}` row of the matrix.
@@ -187,7 +187,7 @@ def lagrange_polynomials(x):
     lagrange_polynomials(4)[0] gives the lagrange polynomials obtained using
     4 LGL points in poly1d form
 
-    lagrange_polynomials(4)[0][2] is :math: `L_2(\\xi)`
+    lagrange_polynomials(4)[0][2] is :math:`L_2(\\xi)`
     lagrange_polynomials(4)[1] gives the coefficients of the above mentioned
     lagrange basis polynomials in a 2D array.
 

dg_maxwell/params.py

@@ -112,7 +112,7 @@
 
 # The wave to be advected is either a sin or a Gaussian wave.
 # This parameter can take values 'sin' or 'gaussian'.
-wave = 'sin'
+wave = 'gaussian'
 
 # Initializing the amplitudes. Change u_init to required initial conditions.
 if (wave=='sin'):

dg_maxwell/tests/test_waveEqn.py

@@ -114,7 +114,7 @@ def change_parameters(LGL, Elements, quad, wave='sin'):
 
 
     # The value of time-step.
-    params.delta_t = params.delta_x / (4 * params.c)
+    params.delta_t = params.delta_x / (4 * abs(params.c))
 
     # Array of timesteps seperated by delta_t.
     params.time    = utils.linspace(0, int(params.total_time / params.delta_t) * params.delta_t,

dg_maxwell/utils.py

@@ -194,6 +194,7 @@ def polyval_1d(polynomials, xi):
     polynomials : af.Array [number_of_polynomials N 1 1]
                  ``number_of_polynomials`` :math:`2D` polynomials of degree
                  :math:`N - 1` of the form
+
                  .. math:: P(x) = a_0x^0 + a_1x^1 + ... \\
                            a_{N - 1}x^{N - 1} + a_Nx^N
     xi      : af.Array [N 1 1 1]

dg_maxwell/wave_equation.py

@@ -49,8 +49,8 @@
 
 def mapping_xi_to_x(x_nodes, xi):
     '''
-    Maps points in :math: `\\xi` space to :math:`x` space using the formula
-    :math:  `x = \\frac{1 - \\xi}{2} x_0 + \\frac{1 + \\xi}{2} x_1`
+    Maps points in :math:`\\xi` space to :math:`x` space using the formula
+    :math:`x = \\frac{1 - \\xi}{2} x_0 + \\frac{1 + \\xi}{2} x_1`
     
     Parameters
     ----------
@@ -59,14 +59,14 @@ def mapping_xi_to_x(x_nodes, xi):
               Element nodes.
     
     xi      : arrayfire.Array [N 1 1 1]
-              Value of :math: `\\xi`coordinate for which the corresponding
-              :math: `x` coordinate is to be found.
+              Value of :math:`\\xi` coordinate for which the corresponding
+              :math:`x` coordinate is to be found.
     
     Returns
     -------
 
     x : arrayfire.Array
-        :math: `x` value in the element corresponding to :math:`\\xi`.
+        :math:`x` value in the element corresponding to :math:`\\xi`.
     '''
 
     N_0 = (1 - xi) / 2
@@ -83,7 +83,7 @@ def mapping_xi_to_x(x_nodes, xi):
 def dx_dxi_numerical(x_nodes, xi):
     '''
 
-    Differential :math: `\\frac{dx}{d \\xi}` calculated by central
+    Differential :math:`\\frac{dx}{d \\xi}` calculated by central
     differential method about xi using the mapping_xi_to_x function.
     
     Parameters
@@ -93,7 +93,7 @@ def dx_dxi_numerical(x_nodes, xi):
               Contains the nodes of elements.
     
     xi      : arrayfire.Array [N_LGL 1 1 1]
-              Values of :math: `\\xi`
+              Values of :math:`\\xi`
     
     Returns
     -------
@@ -115,7 +115,7 @@ def dx_dxi_analytical(x_nodes, xi):
     '''
 
     The analytical result for :math:`\\frac{dx}{d \\xi}` for a 1D element is
-    :math: `\\frac{x_1 - x_0}{2}`
+    :math:`\\frac{x_1 - x_0}{2}`
     
     Parameters
     ----------
@@ -124,14 +124,14 @@ def dx_dxi_analytical(x_nodes, xi):
               Contains the nodes of elements.
  
     xi      : arrayfire.Array [N_LGL 1 1 1]
-              Values of :math: `\\xi`.
+              Values of :math:`\\xi`.
 
     Returns
     -------
 
     analytical_dx_dxi : arrayfire.Array [N_Elements 1 1 1]
-                        The analytical solution of :math:
-                        `\\frac{dx}{d\\xi}` for an element.
+                        The analytical solution of
+                        :math:`\\frac{dx}{d\\xi}` for an element.
     
     '''
     analytical_dx_dxi = (x_nodes[1] - x_nodes[0]) / 2
@@ -143,10 +143,10 @@ def A_matrix():
     '''
 
     Calculates A matrix whose elements :math:`A_{p i}` are given by
-    :math: `A_{p i} &= \\int^1_{-1} L_p(\\xi)L_i(\\xi) \\frac{dx}{d\\xi}`
+    :math:`A_{pi} = \\int^1_{-1} L_p(\\xi)L_i(\\xi) \\frac{dx}{d\\xi}`
 
     The integrals are computed using the integrate() function.
-    Since elements are taken to be of equal size, :math: `\\frac {dx}{dxi}
+    Since elements are taken to be of equal size, :math:`\\frac {dx}{d\\xi}`
     is same everywhere
     
     Returns
@@ -333,8 +333,8 @@ def lax_friedrichs_flux(u_n):
     '''
     Calculates the lax-friedrichs_flux :math:`f_i` using.
 
-    :math:`f_i = \\frac{F(u^{i + 1}_0) + F(u^i_{N_{LGL} - 1})}{2} - \\frac
-                {\Delta x}{2\Delta t} (u^{i + 1}_0 - u^i_{N_{LGL} - 1})`
+    .. math:: f_i = \\frac{F(u^{i + 1}_0) + F(u^i_{N_{LGL} - 1})}{2} - \\
+              \\frac{\Delta x}{2\Delta t} (u^{i + 1}_0 - u^i_{N_{LGL} - 1})
 
     The algorithm used is explained in this `document`_
 
@@ -366,7 +366,72 @@ def lax_friedrichs_flux(u_n):
                         - params.c_lax * (u_iplus1_0 - u_i_N_LGL) / 2
 
 
-    return boundary_flux 
+    return boundary_flux
+
+
+
+def upwind_flux(u_n):
+    '''
+    Finds the upwind flux of all the element edges present inside a domain.
+
+    Parameters
+    ----------
+    u_n : arrayfire.Array [N_LGL N_Elements M 1]
+          Amplitude of the wave at the mapped LGL nodes of each element.
+          This code will work for :math:`M` multiple ``u_n``.
+    
+    Returns
+    -------
+    flux_x : arrayfire.Array [1 N_Elements 1 1]
+             Contains the value of the flux at the boundary elements.
+             Periodic boundary conditions are used.
+    '''
+    right_state    = af.shift(u_n[0, :], 0, -1)
+    left_state     = u_n[-1, :]
+
+    if params.c > 0:
+        return flux_x(left_state)
+
+    if params.c == 0:
+        return flux_x((left_state + right_state) / 2)
+
+    if params.c < 0:
+        return flux_x(right_state)
+
+    return
+
+
+
+def upwind_flux_maxwell_eq(u_n):
+    '''
+    Finds the upwind flux of all the element edges present inside a domain
+    for Mode :math:`1` of Maxwell's equations.
+    Please refer to `maxwell_equation_pbc_1d.pdf`_ for the derivation
+    of the modes.
+
+    .. _maxwell_equation_pbc_1d.pdf: `https://goo.gl/8FS4mv`
+
+    Parameters
+    ----------
+    u_n : arrayfire.Array [N_LGL N_Elements M 1]
+          Amplitude of the wave at the mapped LGL nodes of each element.
+          This code will work for :math:`M` multiple ``u_n``.
+    
+    Returns
+    -------
+    flux_x : arrayfire.Array [1 N_Elements 1 1]
+             Contains the value of the flux at the boundary elements.
+             Periodic boundary conditions are used.
+    '''
+
+    right_state    = af.shift(u_n[0, :], 0, -1)
+
+    flux = right_state.copy()
+
+    flux[:, :, 0] = -right_state[:, :, 1]
+    flux[:, :, 1] = -right_state[:, :, 0]
+
+    return flux
 
 
 
@@ -405,7 +470,15 @@ def surface_term(u_n):
                         d0 = 1, d1 = 1, d2 = shape_u_n[2])
     L_p_1        = af.tile(params.lagrange_basis_value[:, -1],
                         d0 = 1, d1 = 1, d2 = shape_u_n[2])
-    f_i          = lax_friedrichs_flux(u_n)
+    #[NOTE]: Uncomment to use lax friedrichs flux
+    #f_i          = lax_friedrichs_flux(u_n)
+
+    #[NOTE]: Uncomment to use upwind flux for uncoupled advection equations
+    f_i          = upwind_flux(u_n)
+
+    #[NOTE]: Uncomment to use upwind flux for Maxwell's equations
+    #f_i          = upwind_flux_maxwell_eq(u_n)
+
     f_iminus1    = af.shift(f_i, 0, 1)
 
     surface_term = utils.matmul_3D(L_p_1, f_i) \
@@ -483,7 +556,7 @@ def time_evolution(u = None):
     '''
     Solves the wave equation
     
-    ..math: u^{t_n + 1} = b(t_n) \\times A
+    .. math:: u^{t_n + 1} = b(t_n) \\times A
     
     iterated over time.shape[0] time steps t_n 
 
@@ -501,7 +574,7 @@ def time_evolution(u = None):
 
     # Creating a folder to store hdf5 files. If it doesn't exist.
     results_directory = 'results/hdf5_%02d' %(int(params.N_LGL))
-            
+
     if not os.path.exists(results_directory):
         os.makedirs(results_directory)
 
@@ -515,7 +588,7 @@ def time_evolution(u = None):
                         d2 = shape_u_n[2])
 
     element_boundaries = af.np_to_af_array(params.np_element_array)
-    
+
     for t_n in trange(0, time.shape[0]):
         if (t_n % 20) == 0:
             h5file = h5py.File('results/hdf5_%02d/dump_timestep_%06d' %(int(params.N_LGL), int(t_n)) + '.hdf5', 'w')
@@ -524,7 +597,7 @@ def time_evolution(u = None):
             dset[:, :] = u[:, :]
 
         # Code for solving 1D Maxwell's Equations
-        ## Storing u at timesteps which are multiples of 100.
+        # Storing u at timesteps which are multiples of 100.
         #if (t_n % 5) == 0:
             #h5file = h5py.File('results/hdf5_%02d/dump_timestep_%06d' \
                 #%(int(params.N_LGL), int(t_n)) + '.hdf5', 'w')
@@ -536,8 +609,6 @@ def time_evolution(u = None):
             #Ez_dset[:, :] = u[:, :, 0]
             #By_dset[:, :] = u[:, :, 1]
 
-
-        # Implementing RK 4 scheme
         u += RK4_timestepping(A_inverse, u, delta_t)