aboutsummaryrefslogtreecommitdiff
path: root/Examples/Modules/embedded_boundary_rotated_cube/analysis_fields_2d.py
blob: cde45e596f51f0df7fe3a8298114a533e6a138ce (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#! /usr/bin/env python

import yt
import os
import sys
from scipy.constants import mu_0, pi, c
import numpy as np

sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
import checksumAPI

# This is a script that analyses the simulation results from
# the script `inputs_3d`. This simulates a TMmnp mode in a PEC cubic resonator.
# The magnetic field in the simulation is given (in theory) by:
# $$ B_y = \mu \cos(k_x x)\cos(k_z z)\cos( \omega_p t)$$
# with
# $$ k_x = \frac{m\pi}{L}$$
# $$ k_y = \frac{n\pi}{L}$$
# $$ k_z = \frac{p\pi}{L}$$

hi = [0.8, 0.8]
lo = [-0.8, -0.8]
ncells = [32, 32]
dx = (hi[0] - lo[0]) / ncells[0]
dz = (hi[1] - lo[1]) / ncells[1]
m = 0
n = 1
Lx = 1.06
Lz = 1.06

# Open the right plot file
filename = sys.argv[1]
ds = yt.load(filename)
data = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
my_grid = ds.index.grids[0]

By_sim = my_grid['By'].squeeze().v

t = ds.current_time.to_value()

theta = np.pi/8

# Compute the analytic solution
By_th = np.zeros(ncells)
for i in range(ncells[0]):
    for j in range(ncells[1]):
        x = i * dx + lo[0]
        z = j * dz + lo[1]
        xr = x*np.cos(-theta) + z*np.sin(-theta)
        zr = -x*np.sin(-theta) + z*np.cos(-theta)

        By_th[i, j] = mu_0 * (np.cos(m * pi / Lx * (xr - Lx / 2)) *
                              np.cos(n * pi / Lz * (zr - Lz / 2)) *
                              np.cos(np.pi / Lx * c * t))*(By_sim[i, j] != 0)

rel_tol_err = 1e-1

# Compute relative l^2 error on By
rel_err_y = np.sqrt(np.sum(np.square(By_sim - By_th)) / np.sum(np.square(By_th)))
assert (rel_err_y < rel_tol_err)

test_name = os.path.split(os.getcwd())[1]

checksumAPI.evaluate_checksum(test_name, filename)