diff --git a/tetrax/core/experimental_setup.py b/tetrax/core/experimental_setup.py
index ce0f02c847ad61b02f9da0cbe202898585ad8160..6636b4a755fef90daae0ea74e7582289ec06a64a 100644
--- a/tetrax/core/experimental_setup.py
+++ b/tetrax/core/experimental_setup.py
@@ -104,15 +104,13 @@ class ExperimentalSetup:
self.sample.mag = minimize_gibbs_free_energy(self.sample, self._Bext, tol_cg=tol, return_last=return_last,
continue_least_with_squares=continue_least_with_squares)
- def eigenmodes(self, kmin=-40e6, kmax=40e6, Nk=81, num_modes=20, k=None, no_dip=False, h0_dip=False, num_cpus=1,
- save_modes=True,
- save_local=False, save_magphase=False, perturbed_dispersion=False, save_stiffness_fields=False):
+ def eigenmodes(self, kmin=-40e6, kmax=40e6, Nk=81, num_modes=20, k=None, no_dip=False, num_cpus=1, save_modes=True,
+ save_local=False, save_magphase=False):
if self.sample.mag is None:
print("Your sample does not have a magnetization field yet.")
else:
dispersion = self.sample.eigensolver(self.sample, self.Bext, self.name, kmin, kmax, Nk, num_modes, k,
- no_dip, h0_dip, num_cpus, save_modes, save_local, save_magphase,
- perturbed_dispersion, save_stiffness_fields)
+ no_dip, num_cpus, save_modes, save_local, save_magphase)
self._modes_available = save_modes
self.dispersion = dispersion
return dispersion
diff --git a/tetrax/core/operators.py b/tetrax/core/operators.py
index 68a265cb7f1d1c3445e3b4289699b580fd719868..a3ddaeb23a01d21827f1cf98d37e07b1a8a3694f 100644
--- a/tetrax/core/operators.py
+++ b/tetrax/core/operators.py
@@ -187,7 +187,6 @@ class CubicAnisotropyLinearOperator(LinearOperator):
self.Msat = sample.Msat
self.nx = sample.nx
self.m0 = sample.mag
- self.sample = sample
self.fac = 2 * self.Kc / (mu_0 * self.Msat ** 2)
# homogenous anisotropy
@@ -198,15 +197,12 @@ class CubicAnisotropyLinearOperator(LinearOperator):
self.C2 = flattened_mesh_vec_tensor_product(self.c2, self.c2)
self.C3 = flattened_mesh_vec_tensor_product(self.c3, self.c3)
- self.sparse_mat = csr_matrix((3 * self.nx, 3 * self.nx))
-
+ self.sparse_mat = self.make_sparse_mat()
self.shape = self.sparse_mat.shape
- def make_sparse_mat_with_current_m0(self):
- # get current m0
- self.m0 = self.sample.mag.T.flatten()
-
+ def make_sparse_mat(self):
# create bare dyads
+
M0 = flattened_mesh_vec_tensor_product(self.m0, self.m0)
# these are spatially dependent!
@@ -220,9 +216,9 @@ class CubicAnisotropyLinearOperator(LinearOperator):
# A factors are pulled into the tensor product since they are spatially dependent
Nc1 = flattened_mesh_vec_tensor_product(A1 * self.c1, self.c1) + 2 * self.C1.dot(M0).dot(self.C2 + self.C3)
Nc2 = flattened_mesh_vec_tensor_product(A2 * self.c2, self.c2) + 2 * self.C1.dot(M0).dot(self.C2 + self.C3)
- Nc3 = flattened_mesh_vec_tensor_product(A3 * self.c3, self.c3) + 2 * self.C1.dot(M0).dot(self.C2 + self.C3)
+ Nc3 = flattened_mesh_vec_tensor_product(A1 * self.c3, self.c3) + 2 * self.C1.dot(M0).dot(self.C2 + self.C3)
- self.sparse_mat = self.fac * (Nc1 + Nc2 + Nc3)
+ return self.fac * (Nc1 + Nc2 + Nc3)
def set_new_param(self, conf, nx):
"""
diff --git a/tetrax/core/sample.py b/tetrax/core/sample.py
index 3d9834cf650b5df76e2e64d4b0f97f8e99e553d6..85192800f800be0fec3799c528d45d5aaf993185 100644
--- a/tetrax/core/sample.py
+++ b/tetrax/core/sample.py
@@ -6,7 +6,7 @@ import numpy as np
from .fem_core.cythoncore import get_matrices
from .operators import ExchangeOperator, DipolarOperator, UniAxialAnisotropyOperator, BulkDMIOperator, \
- InterfacialDMIOperator, CubicAnisotropyLinearOperator
+ InterfacialDMIOperator
from ..experiments.eigen import calculate_normal_modes, not_implemented_eigensolver
from ..helpers.io import check_mesh_shape
from ..helpers.io import get_mesh_dimension
@@ -158,7 +158,6 @@ class FerromagnetMixin:
self.N_uni = UniAxialAnisotropyOperator(self)
self.N_DMI = BulkDMIOperator(self)
self.N_iDMI = InterfacialDMIOperator(self)
- self.N_cub = CubicAnisotropyLinearOperator(self)
def show(self, show_node_labels=False, comp="vec", scale=5, show_mag=True):
if self.mesh == None:
diff --git a/tetrax/experiments/eigen.py b/tetrax/experiments/eigen.py
index a5f9a941a1046c37b9b5e0333daaff7a80adc923..a65342faaf5bcb138d3f8f52ee64b9a5e1d0fe4a 100644
--- a/tetrax/experiments/eigen.py
+++ b/tetrax/experiments/eigen.py
@@ -7,17 +7,15 @@ import numpy as np
import pandas as pd
import tqdm
from scipy.constants import mu_0
-from scipy.sparse import csc_matrix, csr_matrix
+from scipy.sparse import csc_matrix
from scipy.sparse.linalg import spilu, eigs
from ..core.fem_core.cythoncore import rotation_matrix_py
from ..core.operators import cross_operator, TotalDynMat, SuperLUInv, InvTotalDynMat
-from ..helpers.math import flattened_mesh_vec_scalar_product, diag_matrix_from_vec, flattened_mesh_vec_scalar_product2d, \
- matrix_elem
+from ..helpers.math import flattened_mesh_vec_scalar_product, diag_matrix_from_vec
-def not_implemented_eigensolver(sample, Bext, exp_name, kmin=-40e6, kmax=40e6, Nk=81, num_modes=20, k=None,
- no_dip=False, num_cpus=1, save_modes=False, save_local=False, save_magphase=False):
+def not_implemented_eigensolver(sample, Bext, exp_name, kmin=-40e6, kmax=40e6, Nk=81, num_modes=20, k=None, no_dip=False,num_cpus=1, save_modes=False,save_local=False, save_magphase=False):
print("This eigensolver will be implemented in a future release.")
@@ -26,22 +24,18 @@ def test_rotation(mag, rotation, nx):
Test of rotation matrix and equilibrium vector match, i.e. rotation rotates mag into the z direction.
"""
test_ez = rotation.dot(mag)
- e_z = np.concatenate([np.zeros(2 * nx), np.ones(nx)])
+ e_z = np.concatenate([np.zeros(2*nx), np.ones(nx)])
if not np.allclose(test_ez, e_z):
raise ValueError("Rotation matrix and equilibrium files don't belong together. "
"Please create a new rotation matrix.")
-def calculate_normal_modes(sample, Bext, exp_name, kmin=-40e6, kmax=40e6, Nk=81, num_modes=20, k=None, no_dip=False,
- h0_dip=False,
- num_cpus=1, save_modes=False, save_local=False, save_magphase=False,
- perturbed_dispersion=False, save_stiffness_fields=False):
+def calculate_normal_modes(sample, Bext, exp_name, kmin=-40e6, kmax=40e6, Nk=81, num_modes=20, k=None, no_dip=False,num_cpus=1, save_modes=False,save_local=False, save_magphase=False):
"""
2D version
"""
- experiment_path = "./{}/{}".format(sample.name, exp_name)
- modes_path = "./{}/{}/mode-profiles".format(sample.name, exp_name)
+ modes_path = "./{}/{}/mode-profiles".format(sample.name,exp_name)
if not os.path.exists(modes_path):
os.makedirs(modes_path)
@@ -51,8 +45,7 @@ def calculate_normal_modes(sample, Bext, exp_name, kmin=-40e6, kmax=40e6, Nk=81,
N_uni = sample.N_uni
N_DMI = sample.N_DMI
N_iDMI = sample.N_iDMI
- N_cub = sample.N_cub
- N_cub.make_sparse_mat_with_current_m0()
+
N_iDMI.set_k(0)
N_DMI.set_k(0)
@@ -65,6 +58,7 @@ def calculate_normal_modes(sample, Bext, exp_name, kmin=-40e6, kmax=40e6, Nk=81,
rotation = rotation_matrix_py(mag)
test_rotation(mag, rotation, sample.nx)
+
nx = sample.nx
hext = Bext.T.flatten() / (mu_0 * sample.Msat)
hexc = -N_exc.dot(mag)
@@ -72,29 +66,28 @@ def calculate_normal_modes(sample, Bext, exp_name, kmin=-40e6, kmax=40e6, Nk=81,
hDMI = -N_DMI.dot(mag)
hiDMI = -N_iDMI.dot(mag)
hani = -N_uni.dot(mag)
- hcub = N_cub.nonlinear_field(mag)
+ #hcub = N_cub.nonlinear_field(mag)
if no_dip:
- heff = hexc + hext + hani + hDMI + hiDMI + hcub
- if h0_dip:
- heff += hdip
+ heff = hexc + hext + hani + hDMI + hdip + hiDMI # ++ hcub
else:
- heff = hexc + hdip + hext + hani + hDMI + hiDMI + hcub
+ heff = hexc + hdip + hext + hani + hDMI + hiDMI #+ hcub
h0 = flattened_mesh_vec_scalar_product(heff, mag)
- h0_avrg = np.sum(h0 * sample.dA) / np.sum(sample.dA)
rotation = rotation[:2 * nx, :]
rotation_T = rotation.transpose()
+
# create k and convert to mesh units
if k == None:
- k_ = np.linspace(kmin, kmax, Nk) * sample.scale
+ k_ = np.linspace(kmin,kmax,Nk) * sample.scale
else:
k_ = [k * sample.scale]
- v0 = np.full(2 * nx, 1 + 1j, dtype=complex)
+ v0 = np.full(2*nx, 1+1j, dtype=complex)
+
D_tot = TotalDynMat(N_exc.sparse_mat, N_dip, 1j * mu0_cross[:2 * nx, :2 * nx].dot(rotation), rotation_T)
@@ -102,11 +95,7 @@ def calculate_normal_modes(sample, Bext, exp_name, kmin=-40e6, kmax=40e6, Nk=81,
bndint_order = 6
- modes_per_k_partial = partial(modes_per_k, N_exc, N_uni, N_cub, N_DMI, N_iDMI, h0, h0_avrg, Lambda, rotation,
- rotation_T,
- bndint_order, D_tot, v0, no_dip, sample.scale, sample.Msat, sample.gamma, sample,
- num_modes, modes_path, save_modes, save_local, save_magphase, perturbed_dispersion,
- save_stiffness_fields)
+ modes_per_k_partial = partial(modes_per_k, N_exc, N_uni, N_DMI, N_iDMI, h0, Lambda, rotation, rotation_T, bndint_order, D_tot, v0, no_dip, sample.scale, sample.Msat, sample.gamma, sample, num_modes, modes_path, save_modes,save_local, save_magphase)
if num_cpus == 0 or num_cpus < -1:
num_cpus = 1
@@ -121,39 +110,28 @@ def calculate_normal_modes(sample, Bext, exp_name, kmin=-40e6, kmax=40e6, Nk=81,
else:
pass
+
with mp.Pool(processes=num_cpus) as p:
res_list = []
with tqdm.tqdm(total=len(k_)) as pbar:
for i, res in enumerate(p.imap_unordered(modes_per_k_partial, k_)):
res_list.append(res)
pbar.update()
- df_freqs = pd.concat(res_list, sort=False).sort_values("k (rad/m)").reset_index(drop=True).fillna("nan")
-
- if perturbed_dispersion:
-
- if save_stiffness_fields:
- df_freqs["h0_arvg"] = h0_avrg * np.ones_like(k_)
-
- df_freqs.to_csv(experiment_path + "/dispersion_perturbedmodes.csv")
-
- else:
- df_freqs.to_csv(experiment_path + "/dispersion.csv")
+ df_freqs = pd.concat(res_list,sort=False).sort_values("k (rad/m)").reset_index(drop=True).fillna("nan")
return df_freqs
-def modes_per_k(N_exc, N_uni, N_cub, N_DMI, N_iDMI, h0, h0_avrg, Lambda, rotation, rotation_T, bndint_order, D_tot, v0,
- no_dip, scale, Msat, gamma, sample, num_modes, modes_path, save_modes, save_local, save_magphase,
- perturbed_dispersion, save_stiffness_fields, k):
+def modes_per_k(N_exc, N_uni, N_DMI, N_iDMI, h0, Lambda, rotation, rotation_T, bndint_order, D_tot, v0, no_dip, scale, Msat, gamma,sample, num_modes, modes_path, save_modes, save_local, save_magphase, k):
N_exc.set_k(k)
N_DMI.set_k(k)
N_iDMI.set_k(k)
- N_nodip_k_mat = N_exc.sparse_mat + N_uni.sparse_mat + N_DMI.sparse_mat + N_iDMI.sparse_mat + N_cub.sparse_mat
- h0_operator = diag_matrix_from_vec(np.tile(h0, 3))
- N_nodip_k_mat += h0_operator
+ N_nodip_k_mat = N_exc.sparse_mat + N_uni.sparse_mat + N_DMI.sparse_mat + N_iDMI.sparse_mat
+ N_nodip_k_mat += diag_matrix_from_vec(np.tile(h0, 3))
dyn_mat_k = 1j * Lambda.dot(rotation).dot(N_nodip_k_mat).dot(rotation_T)
+
dyn_mat_k_csc = csc_matrix(dyn_mat_k)
if no_dip:
@@ -179,7 +157,8 @@ def modes_per_k(N_exc, N_uni, N_cub, N_DMI, N_iDMI, h0, h0_avrg, Lambda, rotatio
error_success = this_error.already_sucess
ef = np.zeros(num_modes)
- eigenvectors = np.zeros((2 * D_tot.nx, num_modes))
+ eigenvectors = np.zeros((2*D_tot.nx, num_modes))
+
ef_org = ef.copy()
@@ -195,6 +174,7 @@ def modes_per_k(N_exc, N_uni, N_cub, N_DMI, N_iDMI, h0, h0_avrg, Lambda, rotatio
k_m = k / scale
+
nx = sample.nx
cells = sample.mesh.cells
xyz = sample.xyz
@@ -210,13 +190,13 @@ def modes_per_k(N_exc, N_uni, N_cub, N_DMI, N_iDMI, h0, h0_avrg, Lambda, rotatio
ev_lab_y = ev_lab[nx:-nx]
ev_lab_z = ev_lab[-nx:]
- save_dict = {"Re(m_x)": ev_lab_x.real,
- "Im(m_x)": ev_lab_x.imag,
- "Re(m_y)": ev_lab_y.real,
- "Im(m_y)": ev_lab_y.imag,
- "Re(m_z)": ev_lab_z.real,
+ save_dict = {"Re(m_x)" : ev_lab_x.real,
+ "Im(m_x)" : ev_lab_x.imag,
+ "Re(m_y)" : ev_lab_y.real,
+ "Im(m_y)" : ev_lab_y.imag,
+ "Re(m_z)" : ev_lab_z.real,
"Im(m_z)": ev_lab_z.imag,
- }
+ }
if save_local:
ev_1 = ev[:nx]
@@ -242,7 +222,7 @@ def modes_per_k(N_exc, N_uni, N_cub, N_DMI, N_iDMI, h0, h0_avrg, Lambda, rotatio
save_dict["Arg(m_y)"] = np.angle(ev_lab_y)
save_dict["Arg(m_z)"] = np.angle(ev_lab_z)
- meshio.write_points_cells(filename="{}/mode_k{}radperum_{:03d}.vtk".format(modes_path, k_m * 1e-6, i),
+ meshio.write_points_cells(filename="{}/mode_k{}radperum_{:03d}.vtk".format(modes_path,k_m * 1e-6, i),
points=xyz,
cells=cells,
point_data=save_dict
@@ -253,70 +233,4 @@ def modes_per_k(N_exc, N_uni, N_cub, N_DMI, N_iDMI, h0, h0_avrg, Lambda, rotatio
df_k = pd.DataFrame(columns=["k (rad/m)"] + ["f" + str(j) + " (GHz)" for j in range(len(ef))])
df_k.loc[0] = [k_m] + np.real(ef).tolist()
-
- if perturbed_dispersion:
- # for details of calcultions here, see Phys. Rev. B 104, 174414
- dA = sample.dA
-
- D_tot.set_k(k, csr_matrix((3 * nx, 3 * nx)))
-
- operator_dict = {
- "exc": 1j * Lambda.dot(rotation).dot(N_exc.sparse_mat).dot(rotation_T),
- "uni": 1j * Lambda.dot(rotation).dot(N_uni.sparse_mat).dot(rotation_T),
- "cub": 1j * Lambda.dot(rotation).dot(N_cub.sparse_mat).dot(rotation_T),
- "bDMI": 1j * Lambda.dot(rotation).dot(N_DMI.sparse_mat).dot(rotation_T),
- "iDMI": 1j * Lambda.dot(rotation).dot(N_iDMI.sparse_mat).dot(rotation_T),
- "dip": D_tot
- }
-
- for i in range(eigenvectors.shape[1]):
- Im_N21_ = []
- Re_N21_ = []
- N11_ = []
- N22_ = []
- mode_local = eigenvectors[:, i]
-
- # calculate orthonormal basis from mode profiles
- s1 = np.concatenate([mode_local[:nx], np.zeros(nx)])
- s2 = np.concatenate([np.zeros(nx), mode_local[nx:]])
- s1 /= np.sqrt(np.sum(flattened_mesh_vec_scalar_product2d(np.conjugate(s1), s1) * dA))
- s2 /= np.sqrt(np.sum(flattened_mesh_vec_scalar_product2d(np.conjugate(s2), s2) * dA))
-
- for name in operator_dict:
-
- # in Phys. Rev. B 104, 174414, this is D
- partial_dynmat = operator_dict[name]
-
- # Calculate matrix elements of diagonal part of partial dynamic matrices
- C11 = matrix_elem(s1, s1, partial_dynmat, dA)
- C22 = matrix_elem(s2, s2, partial_dynmat, dA)
- C12 = matrix_elem(s1, s2, partial_dynmat, dA)
- C21 = matrix_elem(s2, s1, partial_dynmat, dA)
-
- # convert to stiffness fields
- Im_N21 = 0.5 * (C11 + C22).real
- Re_N21 = (1j * 0.5 * (C11 - C22)).real
- N11 = (C21).real
- N22 = (C12).real
-
- Im_N21_.append(Im_N21)
- Re_N21_.append(Re_N21)
- N11_.append(N11)
- N22_.append(N22)
-
- if save_stiffness_fields:
- pd.concat([df_k, pd.DataFrame({
- f"Im(N21_{name})_{i}": [Im_N21],
- f"Re(N21_{name})_{i}": [Re_N21],
- f"N11_{name}_{i}": [N11],
- f"N22_{name}_{i}": [N22],
- })], axis=1)
-
- # print(Re_N12)
- omega_i = np.sum(Im_N21_) + np.sqrt(
- (np.sum(N11_) + h0_avrg) * (np.sum(N22_) + h0_avrg) - np.sum(Re_N21_) ** 2)
-
- f_k = omega_i.real * mu_0 * sample.Msat * sample.gamma / 2 / np.pi * 1e-9
- df_k[f"f_pert_{i} (GHz)"] = f_k
-
- return df_k
+ return df_k
\ No newline at end of file
diff --git a/tetrax/experiments/relax.py b/tetrax/experiments/relax.py
index a6d00022cb69a609e2756c8f1db15416c8eea16f..f2cce0a280d6465a2b9ba3a1fee197050fdb95a6 100644
--- a/tetrax/experiments/relax.py
+++ b/tetrax/experiments/relax.py
@@ -28,7 +28,7 @@ def energy_per_length_spherical(mag_spherical, nx, dA, h_ext, N, N_cub):
# print(mag.shape)
# print((h_ext - 0.5 * N.dot(mag)).shape)
# energy_per_volume = -tetramagvec_scalar_product(mag, h_ext + 0.5*N_cub.nonlinear_field(mag) - 0.5 * N.dot(mag)).real
- energy_per_volume = -flattened_mesh_vec_scalar_product(mag, h_ext - 0.5 * N.dot(mag) + 0.5*N_cub.nonlinear_field(mag)).real
+ energy_per_volume = -flattened_mesh_vec_scalar_product(mag, h_ext - 0.5 * N.dot(mag)).real
return (energy_per_volume * dA).sum()
@@ -52,7 +52,7 @@ def jac_spherical(mag_spherical, nx, dA, h_ext, N, N_cub):
m_z = np.cos(theta)
mag = np.concatenate((m_x, m_y, m_z))
- h_eff = h_ext - N.dot(mag).real + N_cub.nonlinear_field(mag)
+ h_eff = h_ext - N.dot(mag).real # + N_cub.nonlinear_field(mag)
# TODO include again!
h_x = h_eff[:nx]
@@ -95,7 +95,7 @@ def minimize_gibbs_free_energy(sample, Bext, tol_cg, return_last=False, continue
sample.N_iDMI.set_k(0)
N_tot = sample.N_dip + sample.N_exc + sample.N_uni + sample.N_DMI + sample.N_iDMI
- N_cub = sample.N_cub
+ N_cub = csr_matrix([])
fac = (sample.Msat * sample.scale) ** 2 * mu_0
# tol_cg = 1e-12
diff --git a/tetrax/helpers/math.py b/tetrax/helpers/math.py
index bd88eec0ce9d13fddb9c4508519d5dec0a916839..31eb071b391744a473b6680b3e7ed005eb8ed38e 100644
--- a/tetrax/helpers/math.py
+++ b/tetrax/helpers/math.py
@@ -8,13 +8,11 @@ is used.
import numpy as np
from scipy.sparse import dia_matrix, bmat
-
def normalize_single_3d_vec(vec):
- norm = np.sqrt(vec.T[0] ** 2 + vec.T[1] ** 2 + vec.T[2] ** 2)
- norm = np.repeat(norm, 3)
- norm = np.reshape(norm, (len(norm) // 3, 3))
- return np.array(vec) / norm
-
+ norm = np.sqrt(vec.T[0]**2 + vec.T[1]**2 + vec.T[2]**2)
+ norm = np.repeat(norm,3)
+ norm = np.reshape(norm,(len(norm)//3,3))
+ return np.array(vec)/norm
def flattened_mesh_vec_abs(vec):
"""Calculates the node-wise magnitude of a flattened mesh vector :math:`\mathbf{v}`.
@@ -136,21 +134,24 @@ def flattened_mesh_vec_tensor_product(vec1, vec2):
tmv2z = vec2[2 * nx:]
# only calculate necessary elements (only 6 distinct) and make diagonal blocks from them
- dxx = dia_matrix((tmv1x * tmv2x, 0), shape=(nx, nx))
- dxy = dia_matrix((tmv1x * tmv2y, 0), shape=(nx, nx))
- dxz = dia_matrix((tmv1x * tmv2z, 0), shape=(nx, nx))
+ dxx = dia_matrix((tmv1x * tmv2x, 0), shape=(nx,nx))
+ dxy = dia_matrix((tmv1x * tmv2y, 0), shape=(nx,nx))
+ dxz = dia_matrix((tmv1x * tmv2z, 0), shape=(nx,nx))
+
+ dyx = dia_matrix((tmv1y * tmv2x, 0), shape=(nx,nx))
+ dyy = dia_matrix((tmv1y * tmv2y, 0), shape=(nx,nx))
+ dyz = dia_matrix((tmv1y * tmv2z, 0), shape=(nx,nx))
- dyx = dia_matrix((tmv1y * tmv2x, 0), shape=(nx, nx))
- dyy = dia_matrix((tmv1y * tmv2y, 0), shape=(nx, nx))
- dyz = dia_matrix((tmv1y * tmv2z, 0), shape=(nx, nx))
+ dzx = dia_matrix((tmv1z * tmv2x, 0), shape=(nx,nx))
+ dzy = dia_matrix((tmv1z * tmv2y, 0), shape=(nx,nx))
+ dzz = dia_matrix((tmv1z * tmv2z, 0), shape=(nx,nx))
- dzx = dia_matrix((tmv1z * tmv2x, 0), shape=(nx, nx))
- dzy = dia_matrix((tmv1z * tmv2y, 0), shape=(nx, nx))
- dzz = dia_matrix((tmv1z * tmv2z, 0), shape=(nx, nx))
# stitch blocks together
- sparse_mat = bmat([[dxx, dxy, dxz], [dyx, dyy, dyz], [dzx, dzy, dzz]], format="csr")
+ sparse_mat = bmat([[dxx, dxy, dxz], [dyx, dyy, dyz], [dzx, dzy, dzz]],format="csr")
+
+
return sparse_mat
@@ -214,8 +215,7 @@ def spherical_angles_to_mesh_vector(theta, phi, as_flattened=False):
vec_y = np.sin(theta) * np.sin(phi)
vec_z = np.cos(theta)
- return np.concatenate((vec_x, vec_y, vec_z)).flatten() if as_flattened else np.concatenate(
- (vec_x, vec_y, vec_z)).flatten().reshape(3, nx).T
+ return np.concatenate((vec_x, vec_y, vec_z)).flatten() if as_flattened else np.concatenate((vec_x, vec_y, vec_z)).flatten().reshape(3, nx).T
def diag_matrix_from_vec(v):
@@ -265,10 +265,10 @@ def h_strip(k, A, L, s):
# sinc is evaluated using np.divide in order to catch divide-by-zero exception
sinc = np.divide(np.sin(k * L / 2), (k * L / 2), out=np.ones_like(k), where=(k != 0))
- return A * sinc * np.exp(-np.abs(k * s))
+ return A*sinc*np.exp(-np.abs(k*s))
-def h_CPW(k, A, L, g, s):
+def h_CPW(k,A,L,g,s):
r"""Fourier transform in `z` direction of one component of the microwave field distribution of a CPW antenna.
Parameters
@@ -307,10 +307,10 @@ def h_CPW(k, A, L, g, s):
Review Research 2, 043203 (2020) <https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.2.043203>`_.
"""
- return 2 * np.sin(k * (g + L) / 2) ** 2 * h_strip(k, A, L, s)
+ return 2 * np.sin(k*(g+L)/2) ** 2 * h_strip(k, A, L, s)
-def h_U(k, A, L, g, s):
+def h_U(k,A,L,g,s):
r"""Fourier transform in `z` direction of one component of the microwave-field distribution of a U-shaped antenna.
Parameters
@@ -349,7 +349,7 @@ def h_U(k, A, L, g, s):
Review Research 2, 043203 (2020) <https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.2.043203>`_.
"""
- return 2 * 1j * np.sin(k * (g + L) / 2) * h_strip(k, A, L, s)
+ return 2 * 1j * np.sin(k*(g+L)/2) * h_strip(k, A, L, s)
def h_hom(k, A):
@@ -379,33 +379,5 @@ def h_hom(k, A):
with :math:`\delta` being the Dirac distribution.
"""
- return A * np.piecewise(k, [k == 0, k != 0], [1, 0])
-
-
-def matrix_elem(i, j, A, dV):
- """
- Calculates the element of an operator `A` in the basis with respect to the basis vectors `i` and `j` as
-
- .. math::
- A_{i,j} = \int\mathrm{d}V \ \mathbf{i}^*\cdot \hat{\mathbf{A}} \cdot \mathbf{j}.
-
- Parameters
- ----------
- i : np.array
- Flattened mesh vector of shape `(3*N,)`.
- j : np.array
- Flattened mesh vector of shape `(3*N,)`.
- A : LinearOperator
- LinearOperator of shape `(3*N, 3*N)`.
- dV : Volume elements of the mesh.
-
- Returns
- -------
- float
- """
- rhs = A.dot(j)
-
- temp = flattened_mesh_vec_scalar_product2d(np.conjugate(i), rhs)
+ return A*np.piecewise(k, [k == 0, k != 0], [1, 0])
- elem = (np.sum(temp * dV))
- return elem