Skip to content
Snippets Groups Projects
Commit cfef25fa authored by Koerber, Lukas (FWIN-C) - 108045's avatar Koerber, Lukas (FWIN-C) - 108045
Browse files

Cleaning up code.

parent c24645e1
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/tetrax/experiments/eigen.py
+ 12
− 126
View file @ cfef25fa
......@@ -55,7 +55,7 @@ def calculate_normal_modes(sample, Bext, exp_name, kmin=-40e6, kmax=40e6, Nk=81,
modes_path = "./{}/{}/mode-profiles".format(sample.name,exp_name)
if not os.path.exists(modes_path):
os.makedirs(modes_path)
#os.mkdir(modes_path)
# exchange operator
N_exc = sample.N_exc
N_dip = sample.N_dip
......@@ -69,33 +69,15 @@ def calculate_normal_modes(sample, Bext, exp_name, kmin=-40e6, kmax=40e6, Nk=81,
N_exc.set_k(0)
N_dip.set_k(0)
tmv_files = ["a"]
#dense_eq = read_tmv_from_project("densemat.dsm", conf)
#dense_eq = np.reshape(dense_eq, (N_dip.nb, N_dip.nb))
#dense_eq_copy = dense_eq.copy()
#N_dip.set_k(0, dense_eq)
#print(float(conf['Msat']))
# exchange length in mesh units
mu0_cross = cross_operator(sample.nx)
#a_n = read_tmv_from_project("area.bin", conf)
#x, y, z = read_coord(conf)
#N_ani = UniAxialAnisotropyOperator(nx, x, y, conf)
#N_cub = CubicAnisotropyLinearOperator(nx, mag, conf)
# rotation = sparsemat_from_files("rotation", conf["project_dir"])
mag = sample.mag.T.flatten()
rotation = rotation_matrix_py(mag)
test_rotation(mag, rotation, sample.nx)
# quit()
nx = sample.nx
hext = Bext.T.flatten() / (mu_0 * sample.Msat) #hom_ext(nx, conf) + hphi_ext(x, y, conf)
hext = Bext.T.flatten() / (mu_0 * sample.Msat)
hexc = -N_exc.dot(mag)
hdip = -N_dip.dot(mag)
hDMI = -N_DMI.dot(mag)
......@@ -104,7 +86,6 @@ def calculate_normal_modes(sample, Bext, exp_name, kmin=-40e6, kmax=40e6, Nk=81,
#hcub = N_cub.nonlinear_field(mag)
if no_dip:
#log.info("Neglecting dipolar interaction.")
heff = hexc + hext + hani + hDMI + hdip + hiDMI # ++ hcub
else:
......@@ -121,53 +102,30 @@ def calculate_normal_modes(sample, Bext, exp_name, kmin=-40e6, kmax=40e6, Nk=81,
k_ = np.linspace(kmin,kmax,Nk) * sample.scale
else:
k_ = [k * sample.scale]
#k_ = [0]
#df_freqs = pd.DataFrame(columns=["k (rad/m)"] + ["f" + str(j) + " (GHz)" for j in range(int(num_modes))])
#df_freqs["k (rad/m)"] = k_
#df_freqs_org = pd.DataFrame(columns=["k (rad/m)"] + ["f" + str(j) + " (GHz)" for j in range(int(num_modes)*2)])
#df_freqs_org["k (rad/m)"] = k_ / sample.scale
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)
#print("we got here")
#quit()
#D_tot = None
#dense_mat_angle = np.diag(read_tmv_from_project(conf["project_name"] + ".sub", conf))
global_idx = 0
Lambda = mu0_cross[:2 * nx, :2 * nx]
bndint_order = 6
#if conf["verbose"]:
#if conf["no_modes"]:
# log.info("not saving any mode profiles")
#else:
# log.info("saving mode profiles")
#log.info("Solving for {} wave vectors and 2*{} modes per k...".format(conf["num_k"], conf["num_modes"]))
#log.verbose("verbose message")
#pbar = tqdm.tqdm(total=len(k_))
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:
#log.warning("You specified an illegal number of CPU cores ({}). Only using 1 instead...".format(conf["num_cpus"]))
num_cpus = 1
if num_cpus == -1:
num_cpus = mp.cpu_count()
if num_cpus > mp.cpu_count():
num_cpus = mp.cpu_count()
#log.warning("You have specified more CPU cores than are available ({}). Using maximum instead...".format(conf["num_cpus"]))
if num_cpus == 1:
#log.info("Using 1 core.")
if mp.cpu_count() > 1:
pass#log.info("Your CPU has {} cores but you are using only 1. Consider using more with the '--cores'. You can use '--cores -1' to utilize all available cores.".format(mp.cpu_count()))
pass
else:
pass#log.info("Using {} cores".format(conf["num_cpus"]))
pass
with mp.Pool(processes=num_cpus) as p:
......@@ -178,33 +136,6 @@ def calculate_normal_modes(sample, Bext, exp_name, kmin=-40e6, kmax=40e6, Nk=81,
pbar.update()
df_freqs = pd.concat(res_list,sort=False).sort_values("k (rad/m)").reset_index(drop=True).fillna("nan")
# create k and convert to mesh units
#if conf["k"] == None:
# if conf["no_dip"]:
# df_freqs.to_csv(str(conf['project_dir'] / "dispersion-eigensolver-no-dip.csv"))
# log.info('Done. Saved dispersion to "dispersion-eigensolver-no-dip.csv"')
# else:
# df_freqs.to_csv(str(conf['project_dir'] / "dispersion-eigensolver.csv"))
# df_freqs_org.to_csv(str(conf['project_dir'] / "dispersion-eigensolver_unsorted.csv"))
# #print("printed_disp")
# log.info('Done. Saved dispersion to "dispersion-eigensolver.csv" and "dispersion-eigensolver_unsorted.csv"')
#else:
# if conf["no_dip"]:
# fname = "frequencies_k{}radperum_nodip.csv".format(conf["k"]*1e-6)
# df_freqs.to_csv(str(conf['project_dir'] / fname))
# log.info('Done. Saved frequencies for k = {k} rad/µm to "{fn}"'.format(
# k=conf["k"]*1e-6,
# fn=fname
# ))
# else:
# fname = "frequencies_k{}radperum.csv".format(conf["k"] * 1e-6)
# df_freqs.to_csv(str(conf['project_dir'] / fname))
# log.info('Done. Saved frequencies for k = {k} rad/µm to "{fn}"'.format(
# k=conf["k"] * 1e-6,
## fn=fname
# ))
return df_freqs
......@@ -217,18 +148,13 @@ def modes_per_k(N_exc, N_uni, N_DMI, N_iDMI, h0, Lambda, rotation, rotation_T, b
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)
#log.verbose("calculating dense matrix on the fly for k*d_tet = {}...".format(k))
#log.verbose("building preconditioner for k*d_tet = {}...".format(k))
#log.verbose("solving for 2*{} modes for k*d_tet = {}...".format(conf["num_modes"],k))
dyn_mat_k_csc = csc_matrix(dyn_mat_k)
if no_dip:
freq, eigenvectors = eigs(dyn_mat_k_csc, v0=v0, which="LM", k=num_modes * 2, tol=1e-3, sigma=0)
ef = freq.real * mu_0 * Msat * gamma / 2 / np.pi * 1e-9
#log.verbose("not including dipolar interaction")
else:
D_tot.set_k(k, N_nodip_k_mat)
try:
......@@ -246,13 +172,10 @@ def modes_per_k(N_exc, N_uni, N_DMI, N_iDMI, h0, Lambda, rotation, rotation_T, b
error_k = this_error.k
error_iter = this_error.after_iter
error_success = this_error.already_sucess
#log.warning("Error inverting dynamic matrix at k = {} : lgmres did not converge (info = {}). "
# "The inverse at this k-value was already calculated successfully {} times. Saving zeros as frequencies".format(error_k, error_iter, error_success))
ef = np.zeros(num_modes)
eigenvectors = np.zeros((2*D_tot.nx, num_modes))
#print("This is the corresponding dense matrix: ")
#print(D_tot.dense)
ef_org = ef.copy()
......@@ -265,32 +188,19 @@ def modes_per_k(N_exc, N_uni, N_DMI, N_iDMI, h0, Lambda, rotation, rotation_T, b
sortindices = np.argsort(ef)
ef = ef[sortindices]
eigenvectors = eigenvectors[:, sortindices]
#log.verbose("positive solutions for k*d_tet = {}:".format(k))
#log.verbose("f (GHz) = " + str(ef))
k_m = k / scale
#if conf["no_dip"]:
# mode_profiles_dir = "mode-profiles-no-dip"
#else:
# mode_profiles_dir = "mode-profiles"
# Path(str(conf['project_dir'] / mode_profiles_dir)).mkdir(parents=True, exist_ok=True)
#fake_i = 0
nx = sample.nx
cells = sample.mesh.cells
xyz = sample.xyz
if save_modes:
#log.verbose("saving mode profiles")
for i in range(eigenvectors.shape[1]):
ev = eigenvectors[:, i]
#if conf["local"]:
# np.real(ev).tofile(str(
# conf['project_dir']) + "/"+ mode_profiles_dir+"/mode_k{}radperum_{:03d}_real_local.tmv".format(k_m * 1e-6,
# fake_i))
# np.imag(ev).tofile(str(
# conf['project_dir']) + "/"+ mode_profiles_dir+"/mode_k{}radperum_{:03d}_imag_local.tmv".format(k_m * 1e-6,
# fake_i))
#log.verbose("Saving mode profiles in local coordinate system, too.")
# rotate eigenvector back to lab system
ev_lab = rotation_T.dot(ev)
ev_lab_x = ev_lab[:nx]
......@@ -329,37 +239,13 @@ def modes_per_k(N_exc, N_uni, N_DMI, N_iDMI, h0, Lambda, rotation, rotation_T, b
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),
points=xyz,
cells=cells,
point_data=save_dict
)
#np.real(ev_lab).tofile(
# str(conf['project_dir'] ) + "/"+ mode_profiles_dir+"/mode_k{}radperum_{:03d}_real.tmv".format(k_m * 1e-6, i))
#np.imag(ev_lab).tofile(
# str(conf['project_dir'] ) + "/"+ mode_profiles_dir+"/mode_k{}radperum_{:03d}_imag.tmv".format(k_m * 1e-6, i))
#if conf["magphase"]:
# #log.verbose("Saving magnitude and phase of mode profiles, too.")
# np.angle(ev_lab).tofile(str(
# conf['project_dir'] ) +"/"+ mode_profiles_dir+"/mode_k{}radperum_{:03d}_phase.tmv".format(
# k_m * 1e-6, fake_i))
# np.abs(ev_lab).tofile(str(
# conf['project_dir'] ) +"/"+ mode_profiles_dir+"/mode_k{}radperum_{:03d}_mag.tmv".format(
# k_m * 1e-6, fake_i))
# if conf["local"]:
# #log.verbose("Saving magnitude and phase in local coordinate system, too")
# np.angle(ev).tofile(str(
# conf['project_dir'] )+"/"+ mode_profiles_dir+"/mode_k{}radperum_{:03d}_phase_local.tmv".format(
# k_m * 1e-6, fake_i))
# np.abs(ev).tofile(str(
# conf['project_dir'] )+"/"+ mode_profiles_dir+"/mode_k{}radperum_{:03d}_mag_local.tmv".format(
# k_m * 1e-6, fake_i))
else:
#log.verbose("not saving any mode profiles for k*d_tet = {}...".format(k))
pass
df_k = pd.DataFrame(columns=["k (rad/m)"] + ["f" + str(j) + " (GHz)" for j in range(len(ef))])
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

A HIFIS Service | Privacy | Imprint | Support | Documentation | Changelog