diff options
Diffstat (limited to 'Examples/Modules')
3 files changed, 195 insertions, 6 deletions
diff --git a/Examples/Modules/embedded_boundary_cube/analysis_fields.py b/Examples/Modules/embedded_boundary_cube/analysis_fields.py index 58dad1bb7..e4b461edd 100755 --- a/Examples/Modules/embedded_boundary_cube/analysis_fields.py +++ b/Examples/Modules/embedded_boundary_cube/analysis_fields.py @@ -32,7 +32,13 @@ Lx = 1 Ly = 1 Lz = 1 h_2 = (m * pi / Lx) ** 2 + (n * pi / Ly) ** 2 + (p * pi / Lz) ** 2 -t = 1.3342563807926085e-08 + +# 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) + +t = ds.current_time.to_value() # Compute the analytic solution Bx_th = np.zeros(ncells) @@ -65,11 +71,6 @@ for i in range(ncells[0]): (-Lz/2 <= z < Lz/2) * np.cos(np.sqrt(2) * np.pi / Lx * c * t)) -# 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) - rel_tol_err = 1e-1 # Compute relative l^2 error on By diff --git a/Examples/Modules/embedded_boundary_rotated_cube/analysis_fields.py b/Examples/Modules/embedded_boundary_rotated_cube/analysis_fields.py new file mode 100755 index 000000000..7765c5d69 --- /dev/null +++ b/Examples/Modules/embedded_boundary_rotated_cube/analysis_fields.py @@ -0,0 +1,116 @@ +#! /usr/bin/env python + +# Copyright 2021 Lorenzo Giacomel +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL + + +import yt +import os, 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 rotated by pi/8. +# The magnetic field in the simulation is given (in theory) by: +# $$ B_x = \frac{-2\mu}{h^2}\, k_x k_z \sin(k_x x)\cos(k_y y)\cos(k_z z)\cos( \omega_p t)$$ +# $$ B_y = \frac{-2\mu}{h^2}\, k_y k_z \cos(k_x x)\sin(k_y y)\cos(k_z z)\cos( \omega_p t)$$ +# $$ B_z = \cos(k_x x)\cos(k_y y)\sin(k_z z)\sin( \omega_p t)$$ +# with +# $$ h^2 = k_x^2 + k_y^2 + k_z^2$$ +# $$ k_x = \frac{m\pi}{L}$$ +# $$ k_y = \frac{n\pi}{L}$$ +# $$ k_z = \frac{p\pi}{L}$$ + +hi = [0.8, 0.8, 0.8] +lo = [-0.8, -0.8, -0.8] +m = 0 +n = 1 +p = 1 +Lx = 1.06 +Ly = 1.06 +Lz = 1.06 +h_2 = (m * pi / Lx) ** 2 + (n * pi / Ly) ** 2 + (p * pi / Lz) ** 2 +theta = np.pi/8 + +# 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) + +t = ds.current_time.to_value() + +rel_tol_err = 1e-2 +my_grid = ds.index.grids[0] + +By_sim = my_grid['raw', 'By_fp'].squeeze().v +Bz_sim = my_grid['raw', 'Bz_fp'].squeeze().v + +ncells = np.array(np.shape(By_sim[:, :, :, 0])) +dx = (hi[0] - lo[0])/ncells[0] +dy = (hi[1] - lo[1])/ncells[1] +dz = (hi[2] - lo[2])/ncells[2] + +# Compute the analytic solution +Bx_th = np.zeros(ncells) +By_th = np.zeros(ncells) +Bz_th = np.zeros(ncells) +for i in range(ncells[0]): + for j in range(ncells[1]): + for k in range(ncells[2]): + x0 = (i+0.5)*dx + lo[0] + y0 = j*dy + lo[1] + z0 = (k+0.5)*dz + lo[2] + + x = x0 + y = y0*np.cos(-theta)-z0*np.sin(-theta) + z = y0*np.sin(-theta)+z0*np.cos(-theta) + By = -2/h_2*mu_0*(n * pi/Ly)*(p * pi/Lz) * (np.cos(m * pi/Lx * (x - Lx/2)) * + np.sin(n * pi/Ly * (y - Ly/2)) * + np.cos(p * pi/Lz * (z - Lz/2)) * + np.cos(np.sqrt(2) * + np.pi / Lx * c * t)) + + Bz = mu_0*(np.cos(m * pi/Lx * (x - Lx/2)) * + np.cos(n * pi/Ly * (y - Ly/2)) * + np.sin(p * pi/Lz * (z - Lz/2)) * + np.cos(np.sqrt(2) * np.pi / Lx * c * t)) + + By_th[i, j, k] = (By*np.cos(theta) - Bz*np.sin(theta))*(By_sim[i, j, k, 0] != 0) + + x0 = (i+0.5)*dx + lo[0] + y0 = (j+0.5)*dy + lo[1] + z0 = k*dz + lo[2] + + x = x0 + y = y0*np.cos(-theta)-z0*np.sin(-theta) + z = y0*np.sin(-theta)+z0*np.cos(-theta) + + By = -2/h_2*mu_0*(n * pi/Ly)*(p * pi/Lz) * (np.cos(m * pi/Lx * (x - Lx/2)) * + np.sin(n * pi/Ly * (y - Ly/2)) * + np.cos(p * pi/Lz * (z - Lz/2)) * + np.cos(np.sqrt(2) * + np.pi / Lx * c * t)) + + Bz = mu_0*(np.cos(m * pi/Lx * (x - Lx/2)) * + np.cos(n * pi/Ly * (y - Ly/2)) * + np.sin(p * pi/Lz * (z - Lz/2)) * + np.cos(np.sqrt(2) * np.pi / Lx * c * t)) + + Bz_th[i, j, k] = (By*np.sin(theta) + Bz*np.cos(theta))*(Bz_sim[i, j, k, 0] != 0) + + +# Compute relative l^2 error on By +rel_err_y = np.sqrt( np.sum(np.square(By_sim[:, :, :, 0] - By_th)) / np.sum(np.square(By_th))) +assert(rel_err_y < rel_tol_err) +# Compute relative l^2 error on Bz +rel_err_z = np.sqrt( np.sum(np.square(Bz_sim[:, :, :, 0] - Bz_th)) / np.sum(np.square(Bz_th))) +assert(rel_err_z < rel_tol_err) + +test_name = os.path.split(os.getcwd())[1] + +checksumAPI.evaluate_checksum(test_name, filename) diff --git a/Examples/Modules/embedded_boundary_rotated_cube/inputs_3d b/Examples/Modules/embedded_boundary_rotated_cube/inputs_3d new file mode 100644 index 000000000..7a6d51cf7 --- /dev/null +++ b/Examples/Modules/embedded_boundary_rotated_cube/inputs_3d @@ -0,0 +1,72 @@ +stop_time = 5.303669113650618e-09 +amr.n_cell = 32 32 32 +amr.max_grid_size = 128 +amr.max_level = 0 + +geometry.coord_sys = 0 +geometry.prob_lo = -0.8 -0.8 -0.8 +geometry.prob_hi = 0.8 0.8 0.8 +warpx.cfl = 1 + +boundary.field_lo = pec pec pec +boundary.field_hi = pec pec pec + +algo.maxwell_solver = ect + +my_constants.xmin = -0.53 +my_constants.ymin = -0.53 +my_constants.zmin = -0.53 +my_constants.xmax = 0.53 +my_constants.ymax = 0.53 +my_constants.zmax = 0.53 +my_constants.pi = 3.141592653589793 +my_constants.theta = pi/8 + +warpx.eb_implicit_function = "max(max(max(x+xmin,-(x+xmax)), max(y*cos(-theta)-z*sin(-theta)+ymin,-(y*cos(-theta)-z*sin(-theta)+ymax))), max(y*sin(-theta)+z*cos(-theta)+zmin,-(y*sin(-theta)+z*cos(-theta)+zmax)))" + +my_constants.m = 0 +my_constants.n = 1 +my_constants.p = 1 +my_constants.Lx = 1.06 +my_constants.Ly = 1.06 +my_constants.Lz = 1.06 +my_constants.x_cent = 0. +my_constants.y_cent = 0. +my_constants.z_cent = 0. +my_constants.h_2 = (m * pi / Lx) ** 2 + (n * pi / Ly) ** 2 + (p * pi / Lz) ** 2 +my_constants.mu_0 = 1.25663706212e-06 + +warpx.B_ext_grid_init_style = parse_B_ext_grid_function +warpx.Bx_external_grid_function(x,y,z) = -2/h_2 * mu_0 * (m * pi / Lx) * (p * pi / Lz) * + sin(m * pi / Lx * (x - Lx / 2 - x_cent)) * + cos(n * pi / Ly * (y*cos(-theta)-z*sin(-theta) - Ly / 2 - y_cent)) * + cos(p * pi / Lz * (y*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent)) + +warpx.By_external_grid_function(x,y,z) = -2/h_2 * mu_0 * (n * pi / Ly) * (p * pi / Lz) * + cos(m * pi / Lx * (x - Lx / 2 - x_cent)) * + sin(n * pi / Ly * (y*cos(-theta)-z*sin(-theta) - Ly / 2 - y_cent)) * + cos(p * pi / Lz * (y*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent)) * + cos(theta) - + mu_0 * + cos(m * pi / Lx * (x - Lx / 2 - x_cent)) * + cos(n * pi / Ly * (y*cos(-theta)-z*sin(-theta) - Ly / 2 - y_cent)) * + sin(p * pi / Lz * (y*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent)) * + sin(theta) + +warpx.Bz_external_grid_function(x,y,z) = mu_0 * + cos(m * pi / Lx * (x - Lx / 2 - x_cent)) * + cos(n * pi / Ly * (y*cos(-theta)-z*sin(-theta) - Ly / 2 - y_cent)) * + sin(p * pi / Lz * (y*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent)) * + cos(theta) - + 2/h_2 * mu_0 * (n * pi / Ly) * (p * pi / Lz) * + cos(m * pi / Lx * (x - Lx / 2 - x_cent)) * + sin(n * pi / Ly * (y*cos(-theta)-z*sin(-theta) - Ly / 2 - y_cent)) * + cos(p * pi / Lz * (y*sin(-theta)+z*cos(-theta) - Lz / 2 - z_cent)) * + sin(theta) + + +diagnostics.diags_names = diag1 +diag1.intervals = 1000 +diag1.diag_type = Full +diag1.plot_raw_fields = 1 +diag1.fields_to_plot = Ex Ey Ez Bx By Bz |