Developmental file for modifying the 1D advection solver to work for multiple wave equations


In [1]:
import os
import sys
sys.path.insert(0, os.path.abspath('../../'))

import numpy as np
from matplotlib import pyplot as plt
import arrayfire as af

from dg_maxwell import params
from dg_maxwell import lagrange
from dg_maxwell import wave_equation as w1d
from dg_maxwell import utils

af.set_backend('opencl')
af.set_device(1)
af.info()

plt.rcParams['figure.figsize']     = 12, 7.5
plt.rcParams['lines.linewidth']    = 1.5
plt.rcParams['font.family']        = 'serif'
plt.rcParams['font.weight']        = 'bold'
plt.rcParams['font.size']          = 20  
plt.rcParams['font.sans-serif']    = 'serif'
plt.rcParams['text.usetex']        = True
plt.rcParams['axes.linewidth']     = 1.5
plt.rcParams['axes.titlesize']     = 'medium'
plt.rcParams['axes.labelsize']     = 'medium'

plt.rcParams['xtick.major.size']   = 8
plt.rcParams['xtick.minor.size']   = 4
plt.rcParams['xtick.major.pad']    = 8
plt.rcParams['xtick.minor.pad']    = 8
plt.rcParams['xtick.color']        = 'k'
plt.rcParams['xtick.labelsize']    = 'medium'
plt.rcParams['xtick.direction']    = 'in'    

plt.rcParams['ytick.major.size']   = 8
plt.rcParams['ytick.minor.size']   = 4
plt.rcParams['ytick.major.pad']    = 8
plt.rcParams['ytick.minor.pad']    = 8
plt.rcParams['ytick.color']        = 'k'
plt.rcParams['ytick.labelsize']    = 'medium'
plt.rcParams['ytick.direction']    = 'in'
plt.rcParams['text.usetex']        = True
plt.rcParams['text.latex.unicode'] = True


/home/ubermensch/.local/anaconda3/lib/python3.6/site-packages/numpy/lib/polynomial.py:1193: FutureWarning: In the future extra properties will not be copied across when constructing one poly1d from another
  other = poly1d(other)
/home/ubermensch/.local/anaconda3/lib/python3.6/site-packages/numpy/lib/polynomial.py:1220: FutureWarning: In the future extra properties will not be copied across when constructing one poly1d from another
  other = poly1d(other)

In [2]:
# 1. Set the initial conditions

E_00 = 1.
E_01 = 1.

B_00 = 0.2
B_01 = 0.5

E_z_init = E_00 * af.sin(2 * np.pi * params.element_LGL) \
         + E_01 * af.cos(2 * np.pi * params.element_LGL)

B_y_init = B_00 * af.sin(2 * np.pi * params.element_LGL) \
         + B_01 * af.cos(2 * np.pi * params.element_LGL)

u_init = af.constant(0., d0 = params.N_LGL, d1 = params.N_Elements, d2 = 2)
u_init[:, :, 0] = E_z_init
u_init[:, :, 1] = B_y_init

In [3]:
element_LGL_flat = af.flat(params.element_LGL)
E_z_init_flat    = af.flat(u_init[:, :, 0])
B_y_init_flat    = af.flat(u_init[:, :, 1])

plt.plot(element_LGL_flat, E_z_init_flat, label = r'$E_z$')
plt.plot(element_LGL_flat, B_y_init_flat, label = r'$B_y$')

plt.title(r'Plot of $E_z(t = 0)$ and $B_y(t = 0)$')
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')

plt.legend(prop={'size': 14})

plt.show()


Prototype implementation of RK4 for multiple-u's


In [20]:
u = u_init[:, :, :]

shape_u_n = utils.shape(u)

In [21]:
A_inverse = af.tile(af.inverse(w1d.A_matrix()),
                    d0 = 1, d1 = 1,
                    d2 = shape_u_n[2])


k1 = utils.matmul_3D(A_inverse, w1d.b_vector(u))
k2 = utils.matmul_3D(A_inverse, w1d.b_vector(u + k1 * params.delta_t / 2))
k3 = utils.matmul_3D(A_inverse, w1d.b_vector(u + k2 * params.delta_t / 2))
k4 = utils.matmul_3D(A_inverse, w1d.b_vector(u + k3 * params.delta_t))

delta_u = params.delta_t * (k1 + 2 * k2 + 2 * k3 + k4) / 6

# print(delta_u)

In [3]:
w1d.time_evolution(u_init)


100%|██████████| 626/626 [00:46<00:00, 13.54it/s]
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-3-e6442e9e8d4d> in <module>()
----> 1 w1d.time_evolution(u_init)

~/workspace/quazar/git/amanabt/DG_Maxwell/dg_maxwell/wave_equation.py in time_evolution(u)
    805     u_analytical = analytical_u_LGL(t_n + 1)
    806 
--> 807     u_diff = af.abs(u - u_analytical)
    808 
    809 

~/.local/anaconda3/lib/python3.6/site-packages/arrayfire/array.py in __sub__(self, other)
    850         Return self - other.
    851         """
--> 852         return _binary_func(self, other, backend.get().af_sub)
    853 
    854     def __isub__(self, other):

~/.local/anaconda3/lib/python3.6/site-packages/arrayfire/array.py in _binary_func(lhs, rhs, c_func)
    173         raise TypeError("Invalid parameter to binary function")
    174 
--> 175     safe_call(c_func(c_pointer(out.arr), lhs.arr, other.arr, _bcast_var.get()))
    176 
    177     return out

~/.local/anaconda3/lib/python3.6/site-packages/arrayfire/util.py in safe_call(af_error)
     77         err_len = c_dim_t(0)
     78         backend.get().af_get_last_error(c_pointer(err_str), c_pointer(err_len))
---> 79         raise RuntimeError(to_str(err_str))
     80 
     81 def get_version():

RuntimeError: In function af::dim4 getOutDims(const af::dim4&, const af::dim4&, bool)
In file src/backend/ArrayInfo.cpp:178
Invalid dimension for argument 1
Expected: ldims == rdims

In [ ]: