Skip to content
GitLab
Explore
Sign in
Primary navigation
Search or go to…
Project
TetraX
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Deploy
Releases
Package Registry
Model registry
Operate
Terraform modules
Monitor
Service Desk
Analyze
Contributor analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Terms and privacy
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
Micromagnetic modeling
TetraX
Commits
e7769169
Commit
e7769169
authored
3 years ago
by
Koerber, Lukas (FWIN-C) - 108045
Browse files
Options
Downloads
Patches
Plain Diff
Implement possibility to set multiple (and inhomogeneous) uniaxial anistropies.
parent
5f1e8f71
No related branches found
Branches containing commit
No related tags found
Tags containing commit
No related merge requests found
Changes
2
Hide whitespace changes
Inline
Side-by-side
Showing
2 changed files
src/tetrax/core/operators.py
+126
-118
126 additions, 118 deletions
src/tetrax/core/operators.py
src/tetrax/core/sample.py
+3
-0
3 additions, 0 deletions
src/tetrax/core/sample.py
with
129 additions
and
118 deletions
src/tetrax/core/operators.py
+
126
−
118
View file @
e7769169
#from tetraeigen2d_lib.IO_helpers import sparsemat_from_project, read_tmv_from_project, read_coord
#
from tetraeigen2d_lib.IO_helpers import sparsemat_from_project, read_tmv_from_project, read_coord
from
scipy.sparse.linalg
import
spsolve
,
LinearOperator
,
gmres
,
lgmres
,
cgs
from
scipy.sparse
import
csr_matrix
,
dia_matrix
,
coo_matrix
,
identity
,
csc_matrix
,
diags
,
bmat
import
numpy
as
np
...
...
@@ -13,18 +13,19 @@ class ExchangeOperator(LinearOperator):
def
__init__
(
self
,
sample
):
self
.
sparse_mat_k0
=
sample
.
xc
#self.sparse_mat = -xc_op_py("{}.matricesandmisc/{}".format(conf["project_name"],conf["project_name"]))
#
self.sparse_mat = -xc_op_py("{}.matricesandmisc/{}".format(conf["project_name"],conf["project_name"]))
self
.
shape
=
self
.
sparse_mat_k0
.
shape
#print(self.shape)
#
print(self.shape)
self
.
dtype
=
np
.
complex
self
.
lex
=
np
.
sqrt
(
2
*
sample
.
Aex
/
(
mu_0
*
sample
.
Msat
**
2
)
)
/
sample
.
scale
self
.
lex
=
np
.
sqrt
(
2
*
sample
.
Aex
/
(
mu_0
*
sample
.
Msat
**
2
))
/
sample
.
scale
self
.
k
=
0
self
.
nx
=
sample
.
nx
self
.
set_k
(
0
)
def
set_k
(
self
,
k
):
self
.
k
=
k
self
.
sparse_mat
=
self
.
lex
**
2
*
self
.
sparse_mat_k0
+
self
.
lex
**
2
*
self
.
k
**
2
*
dia_matrix
((
np
.
ones
(
3
*
self
.
nx
),
0
),
shape
=
self
.
shape
)
self
.
sparse_mat
=
self
.
lex
**
2
*
self
.
sparse_mat_k0
+
self
.
lex
**
2
*
self
.
k
**
2
*
dia_matrix
(
(
np
.
ones
(
3
*
self
.
nx
),
0
),
shape
=
self
.
shape
)
def
_matvec
(
self
,
vec
):
return
self
.
sparse_mat
.
dot
(
vec
)
...
...
@@ -33,23 +34,24 @@ class ExchangeOperator(LinearOperator):
class
BulkDMIOperator
(
LinearOperator
):
def
__init__
(
self
,
sample
):
self
.
grad_x
=
sample
.
grad_x
self
.
grad_y
=
sample
.
grad_y
self
.
nx
=
sample
.
nx
self
.
fac
=
2
*
sample
.
Dbulk
*
(
sample
.
Msat
**
2
*
mu_0
*
sample
.
scale
)
self
.
fac
=
2
*
sample
.
Dbulk
*
(
sample
.
Msat
**
2
*
mu_0
*
sample
.
scale
)
self
.
k
=
0
# construct Sigma matrix (for k-dep part)
Sigma_mat_x
=
np
.
hstack
((
np
.
zeros
((
self
.
nx
,
self
.
nx
)),
-
1j
*
np
.
diag
(
np
.
ones
(
self
.
nx
)),
np
.
zeros
((
self
.
nx
,
self
.
nx
))))
Sigma_mat_y
=
np
.
hstack
((
1j
*
np
.
diag
(
np
.
ones
(
self
.
nx
)),
np
.
zeros
((
self
.
nx
,
self
.
nx
)),
np
.
zeros
((
self
.
nx
,
self
.
nx
))))
Sigma_mat_z
=
np
.
hstack
((
np
.
zeros
((
self
.
nx
,
self
.
nx
)),
np
.
zeros
((
self
.
nx
,
self
.
nx
)),
np
.
zeros
((
self
.
nx
,
self
.
nx
))))
Sigma_mat_x
=
np
.
hstack
(
(
np
.
zeros
((
self
.
nx
,
self
.
nx
)),
-
1j
*
np
.
diag
(
np
.
ones
(
self
.
nx
)),
np
.
zeros
((
self
.
nx
,
self
.
nx
))))
Sigma_mat_y
=
np
.
hstack
(
(
1j
*
np
.
diag
(
np
.
ones
(
self
.
nx
)),
np
.
zeros
((
self
.
nx
,
self
.
nx
)),
np
.
zeros
((
self
.
nx
,
self
.
nx
))))
Sigma_mat_z
=
np
.
hstack
(
(
np
.
zeros
((
self
.
nx
,
self
.
nx
)),
np
.
zeros
((
self
.
nx
,
self
.
nx
)),
np
.
zeros
((
self
.
nx
,
self
.
nx
))))
Sigma_mat
=
np
.
vstack
((
Sigma_mat_x
,
Sigma_mat_y
,
Sigma_mat_z
))
self
.
Sigma
=
csr_matrix
(
Sigma_mat
)
# construct Sigma matrix (for k-dep part)
rot_mat_x
=
np
.
hstack
((
np
.
zeros
((
self
.
nx
,
self
.
nx
)),
np
.
zeros
((
self
.
nx
,
self
.
nx
)),
self
.
grad_y
.
todense
()))
rot_mat_y
=
np
.
hstack
((
np
.
zeros
((
self
.
nx
,
self
.
nx
)),
np
.
zeros
((
self
.
nx
,
self
.
nx
)),
-
self
.
grad_x
.
todense
()))
...
...
@@ -59,11 +61,11 @@ class BulkDMIOperator(LinearOperator):
self
.
rot
=
csr_matrix
(
rot_mat
)
self
.
sparse_mat
=
self
.
fac
*
(
self
.
rot
+
self
.
k
*
self
.
Sigma
)
self
.
sparse_mat
=
self
.
fac
*
(
self
.
rot
+
self
.
k
*
self
.
Sigma
)
self
.
shape
=
self
.
sparse_mat
.
shape
def
set_k
(
self
,
k
):
def
set_k
(
self
,
k
):
self
.
k
=
k
self
.
sparse_mat
=
self
.
fac
*
(
self
.
rot
+
self
.
k
*
self
.
Sigma
)
...
...
@@ -74,49 +76,42 @@ class BulkDMIOperator(LinearOperator):
class
InterfacialDMIOperator
(
LinearOperator
):
def
__init__
(
self
,
sample
):
self
.
grad_x
=
sample
.
grad_x
self
.
grad_y
=
sample
.
grad_y
self
.
nx
=
sample
.
nx
# for now symmetry breaking direction is in y.
self
.
d_vec
=
np
.
repeat
(
np
.
array
([
0
,
1
,
0
]),
self
.
nx
,
axis
=
0
).
flatten
()
self
.
fac
=
2
*
sample
.
Didmi
*
(
sample
.
Msat
**
2
*
mu_0
*
sample
.
scale
)
self
.
nx
,
axis
=
0
).
flatten
()
self
.
fac
=
2
*
sample
.
Didmi
*
(
sample
.
Msat
**
2
*
mu_0
*
sample
.
scale
)
self
.
k
=
0
# construct Pi matrix (for k-dep part)
d_vec_x
=
self
.
d_vec
[:
self
.
nx
]
d_vec_x_diag
=
dia_matrix
((
d_vec_x
,
0
),
shape
=
(
self
.
nx
,
self
.
nx
))
d_vec_x_diag
=
dia_matrix
((
d_vec_x
,
0
),
shape
=
(
self
.
nx
,
self
.
nx
))
d_vec_y
=
self
.
d_vec
[
self
.
nx
:
-
self
.
nx
]
d_vec_y_diag
=
dia_matrix
((
d_vec_y
,
0
),
shape
=
(
self
.
nx
,
self
.
nx
))
zeros_block
=
csr_matrix
((
self
.
nx
,
self
.
nx
))
self
.
Pi_mat
=
bmat
([[
zeros_block
,
zeros_block
,
1j
*
d_vec_x_diag
],
[
zeros_block
,
zeros_block
,
1j
*
d_vec_y_diag
],
[
-
1j
*
d_vec_x_diag
,
-
1j
*
d_vec_y_diag
,
zeros_block
]])
zeros_block
=
csr_matrix
((
self
.
nx
,
self
.
nx
))
self
.
Pi_mat
=
bmat
([[
zeros_block
,
zeros_block
,
1j
*
d_vec_x_diag
],
[
zeros_block
,
zeros_block
,
1j
*
d_vec_y_diag
],
[
-
1j
*
d_vec_x_diag
,
-
1j
*
d_vec_y_diag
,
zeros_block
]])
self
.
diff_mat
=
bmat
([[
d_vec_x_diag
.
dot
(
self
.
grad_x
),
d_vec_x_diag
.
dot
(
self
.
grad_y
),
zeros_block
],
[
d_vec_y_diag
.
dot
(
self
.
grad_x
),
d_vec_y_diag
.
dot
(
self
.
grad_y
),
zeros_block
],
[
zeros_block
,
zeros_block
,
zeros_block
]])
-
\
[
zeros_block
,
zeros_block
,
zeros_block
]])
-
\
bmat
([[
self
.
grad_x
.
dot
(
d_vec_x_diag
),
self
.
grad_x
.
dot
(
d_vec_y_diag
),
zeros_block
],
[
self
.
grad_y
.
dot
(
d_vec_x_diag
),
self
.
grad_y
.
dot
(
d_vec_y_diag
),
zeros_block
],
[
zeros_block
,
zeros_block
,
zeros_block
]])
self
.
sparse_mat
=
self
.
fac
*
(
self
.
diff_mat
+
self
.
k
*
self
.
Pi_mat
)
self
.
sparse_mat
=
self
.
fac
*
(
self
.
diff_mat
+
self
.
k
*
self
.
Pi_mat
)
self
.
shape
=
self
.
sparse_mat
.
shape
def
set_k
(
self
,
k
):
def
set_k
(
self
,
k
):
self
.
k
=
k
self
.
sparse_mat
=
self
.
fac
*
(
self
.
diff_mat
+
self
.
k
*
self
.
Pi_mat
)
self
.
sparse_mat
=
self
.
fac
*
(
self
.
diff_mat
+
self
.
k
*
self
.
Pi_mat
)
def
_matvec
(
self
,
x
):
return
self
.
sparse_mat
.
dot
(
x
)
...
...
@@ -124,34 +119,51 @@ class InterfacialDMIOperator(LinearOperator):
class
UniAxialAnisotropyOperator
(
LinearOperator
):
def
__init__
(
self
,
sample
):
def
__init__
(
self
,
sample
)
->
None
:
self
.
fac
=
-
2
*
sample
.
Ku1
/
(
mu_0
*
sample
.
Msat
**
2
)
self
.
k_theta
=
sample
.
k_theta
self
.
k_phi
=
sample
.
k_phi
self
.
Msat
=
sample
.
Msat
self
.
xyz_shape
=
sample
.
xyz
.
shape
self
.
Ku1_list
=
[
sample
.
Ku1
]
if
isinstance
(
sample
.
Ku1
,
(
float
,
int
,
np
.
ndarray
))
else
sample
.
Ku1
self
.
e_u_list
=
[
sample
.
e_u
]
if
isinstance
(
sample
.
e_u
,
np
.
ndarray
)
else
sample
.
e_u
self
.
nx
=
sample
.
nx
self
.
sparse_mat
=
self
.
make_sparse_mat
()
self
.
shape
=
self
.
sparse_mat
.
shape
self
.
dtype
=
np
.
complex128
if
len
(
self
.
Ku1_list
)
!=
len
(
self
.
e_u_list
):
print
(
f
"
Could not initialize UniAxialAnisotropyOperator. Not the same number of Ku1 (
{
len
(
self
.
Ku1_list
)
}
) and e_u (
{
len
(
self
.
e_u_list
)
}
) specified.
"
)
else
:
self
.
sparse_mat
=
self
.
make_sparse_mat
()
self
.
shape
=
self
.
sparse_mat
.
shape
self
.
dtype
=
np
.
complex128
def
make_sparse_mat
(
self
):
sparse_mat
=
csr_matrix
((
3
*
self
.
nx
,
3
*
self
.
nx
))
# iterate over all anisotropies.
for
i
,
Ku1
in
enumerate
(
self
.
Ku1_list
):
# make a flattened mesh vector Ku1*ones(3*nx) from Ku1 if it is a single number,
# tile (Ku1, Ku1, Ku1) to mesh vector if inhomogeneous (array).
Ku1_array
=
np
.
tile
(
Ku1
,
3
)
if
isinstance
(
Ku1
,
np
.
ndarray
)
else
Ku1
*
np
.
ones
(
3
*
self
.
nx
)
theta
=
self
.
k_theta
/
180
*
np
.
pi
phi
=
self
.
k_phi
/
180
*
np
.
pi
k_vec
=
np
.
array
([
np
.
sin
(
theta
)
*
np
.
cos
(
phi
),
np
.
sin
(
theta
)
*
np
.
sin
(
phi
),
np
.
cos
(
theta
)])
e_u
=
self
.
e_u_list
[
i
]
k_vec_
=
np
.
repeat
(
k_vec
,
self
.
nx
,
axis
=
0
).
flatten
()
# check if e_u is single vector or inhomogeneous
if
np
.
shape
(
e_u
)
==
self
.
xyz_shape
:
# only normalize and convert to flattened mesh vector
e_u_normalized
=
normalize_single_3d_vec
(
e_u
).
T
.
flatten
()
sparse_mat
=
flattened_mesh_vec_tensor_product
(
k_vec_
,
k_vec_
)
elif
np
.
shape
(
e_u
)
==
(
3
,):
return
self
.
fac
*
sparse_mat
# extent vector to mesh vector, normalize and flatten
# TODO this could be faster by switching the order of things
e_u_normalized
=
normalize_single_3d_vec
(
np
.
array
([
e_u
for
i
in
range
(
self
.
nx
)])).
T
.
flatten
()
sparse_mat
+=
flattened_mesh_vec_tensor_product
(
-
2
*
Ku1_array
/
(
mu_0
*
self
.
Msat
**
2
)
*
e_u_normalized
,
e_u_normalized
)
return
sparse_mat
def
_matvec
(
self
,
vec
):
return
self
.
sparse_mat
.
dot
(
vec
)
...
...
@@ -159,13 +171,12 @@ class UniAxialAnisotropyOperator(LinearOperator):
class
CubicAnisotropyLinearOperator
(
LinearOperator
):
def
__init__
(
self
,
sample
):
self
.
Kc
=
sample
.
Kc1
self
.
Msat
=
sample
.
Msat
self
.
nx
=
sample
.
nx
self
.
m0
=
sample
.
mag
self
.
fac
=
2
*
self
.
Kc
/
(
mu_0
*
self
.
Msat
**
2
)
self
.
fac
=
2
*
self
.
Kc
/
(
mu_0
*
self
.
Msat
**
2
)
# homogenous anisotropy
self
.
c1
=
np
.
repeat
(
np
.
array
(
sample
.
v1Kc
),
self
.
nx
,
axis
=
0
).
flatten
()
...
...
@@ -175,28 +186,28 @@ 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
=
self
.
make_sparse_mat
()
self
.
shape
=
self
.
sparse_mat
.
shape
def
make_sparse_mat
(
self
):
# create bare dyads
M0
=
flattened_mesh_vec_tensor_product
(
self
.
m0
,
self
.
m0
)
# these are spatially dependent!
A1
=
np
.
tile
((
flattened_mesh_vec_scalar_product
(
self
.
m0
,
self
.
c2
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
self
.
m0
,
self
.
c3
))
**
2
,
3
).
flatten
()
A2
=
np
.
tile
((
flattened_mesh_vec_scalar_product
(
self
.
m0
,
self
.
c1
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
self
.
m0
,
self
.
c3
))
**
2
,
3
).
flatten
()
A3
=
np
.
tile
((
flattened_mesh_vec_scalar_product
(
self
.
m0
,
self
.
c1
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
self
.
m0
,
self
.
c2
))
**
2
,
3
).
flatten
()
A1
=
np
.
tile
((
flattened_mesh_vec_scalar_product
(
self
.
m0
,
self
.
c2
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
self
.
m0
,
self
.
c3
))
**
2
,
3
).
flatten
()
A2
=
np
.
tile
((
flattened_mesh_vec_scalar_product
(
self
.
m0
,
self
.
c1
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
self
.
m0
,
self
.
c3
))
**
2
,
3
).
flatten
()
A3
=
np
.
tile
((
flattened_mesh_vec_scalar_product
(
self
.
m0
,
self
.
c1
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
self
.
m0
,
self
.
c2
))
**
2
,
3
).
flatten
()
# 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
(
A1
*
self
.
c3
,
self
.
c3
)
+
2
*
self
.
C1
.
dot
(
M0
).
dot
(
self
.
C2
+
self
.
C3
)
return
self
.
fac
*
(
Nc1
+
Nc2
+
Nc3
)
return
self
.
fac
*
(
Nc1
+
Nc2
+
Nc3
)
def
set_new_param
(
self
,
conf
,
nx
):
"""
...
...
@@ -205,16 +216,19 @@ class CubicAnisotropyLinearOperator(LinearOperator):
raise
NotImplementedError
def
_matvec
(
self
,
vec
):
#print(self.sparse_mat)
#print(vec)
#
print(self.sparse_mat)
#
print(vec)
return
self
.
sparse_mat
.
dot
(
vec
)
def
nonlinear_field_old
(
self
,
vec
):
def
nonlinear_field_old
(
self
,
vec
):
"""
Full nonlinear anistropy. Returns the cubic anisotropy field of a given vector.
This is not a linear operator.
"""
A1
=
np
.
tile
((
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c2
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c3
))
**
2
,
3
).
flatten
()
A2
=
np
.
tile
((
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c1
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c3
))
**
2
,
3
).
flatten
()
A3
=
np
.
tile
((
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c1
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c2
))
**
2
,
3
).
flatten
()
A1
=
np
.
tile
((
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c2
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c3
))
**
2
,
3
).
flatten
()
A2
=
np
.
tile
((
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c1
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c3
))
**
2
,
3
).
flatten
()
A3
=
np
.
tile
((
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c1
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c2
))
**
2
,
3
).
flatten
()
hc1
=
flattened_mesh_vec_tensor_product
(
A1
*
self
.
c1
,
self
.
c1
).
dot
(
vec
)
hc2
=
flattened_mesh_vec_tensor_product
(
A2
*
self
.
c2
,
self
.
c2
).
dot
(
vec
)
...
...
@@ -227,45 +241,44 @@ class CubicAnisotropyLinearOperator(LinearOperator):
"""
Full nonlinear anistropy. Returns the cubic anisotropy field of a given vector.
This is not a linear operator.
"""
A1
=
np
.
tile
(
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c2
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c3
))
**
2
,
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c2
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c3
))
**
2
,
3
).
flatten
()
A2
=
np
.
tile
(
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c1
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c3
))
**
2
,
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c1
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c3
))
**
2
,
3
).
flatten
()
A3
=
np
.
tile
(
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c1
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c2
))
**
2
,
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c1
))
**
2
+
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c2
))
**
2
,
3
).
flatten
()
hc1
=
A1
*
np
.
tile
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c1
),
3
).
flatten
()
*
self
.
c1
hc2
=
A2
*
np
.
tile
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c2
),
3
).
flatten
()
*
self
.
c2
hc3
=
A3
*
np
.
tile
(
flattened_mesh_vec_scalar_product
(
vec
,
self
.
c3
),
3
).
flatten
()
*
self
.
c3
hc
=
-
self
.
fac
*
(
hc1
+
hc2
+
hc3
)
return
hc
def
N_dip
(
m
,
k
,
ksquared_mat_a_n
,
ksquared_mat_a_nb
,
nx
,
div_x
,
div_y
,
poiss
,
boundary_nodes
,
dense
,
laplace
,
grad_x
,
grad_y
,
a_n
):
div_x
,
div_y
,
poiss
,
boundary_nodes
,
dense
,
laplace
,
grad_x
,
grad_y
,
a_n
):
# calculate charges
m_x
=
m
[:
nx
]
m_y
=
m
[
nx
:
2
*
nx
]
m_z
=
m
[
2
*
nx
:]
m_y
=
m
[
nx
:
2
*
nx
]
m_z
=
m
[
2
*
nx
:]
rhs
=
div_x
.
dot
(
m_x
)
+
div_y
.
dot
(
m_y
)
+
1j
*
k
*
m_z
*
a_n
rhs
=
div_x
.
dot
(
m_x
)
+
div_y
.
dot
(
m_y
)
+
1j
*
k
*
m_z
*
a_n
# (nable^2 - k^2) * psi1_k = rho_k
#op = poiss - ksquared_mat_a_n*3
#if 0.015 < k < 0.03:
#
op = poiss - ksquared_mat_a_n*3
#
if 0.015 < k < 0.03:
# print(k, det(op.todense()))
#rho_k2 = -(poiss - ksquared_mat_a_n).dot(rho_k)
#
rho_k2 = -(poiss - ksquared_mat_a_n).dot(rho_k)
psi1
=
spsolve
(
poiss
-
ksquared_mat_a_n
,
rhs
)
# get boundary conditions for psi2
sigma
=
np
.
zeros
(
nx
)
+
0j
...
...
@@ -278,16 +291,15 @@ def N_dip(m, k, ksquared_mat_a_n, ksquared_mat_a_nb, nx,
# solve for psi2 potential
psi2
=
spsolve
(
laplace
-
ksquared_mat_a_n
,
-
sigma2
)
# add up to make full potential
psi
=
psi1
+
psi2
psi
=
psi1
+
psi2
# calculate (negative) lateral dipolar field
Nm_x
=
grad_x
.
dot
(
psi
)
Nm_y
=
grad_y
.
dot
(
psi
)
# (negative) longitudinal field
Nm_z
=
1j
*
k
*
psi
Nm_z
=
1j
*
k
*
psi
return
np
.
array
([
Nm_x
,
Nm_y
,
Nm_z
]).
flatten
()
...
...
@@ -308,14 +320,14 @@ class DipolarOperator(LinearOperator):
self
.
ksquared_mat
=
None
self
.
nx
=
sample
.
nx
self
.
nb
=
sample
.
nb
self
.
shape
=
3
*
self
.
nx
,
3
*
self
.
nx
self
.
shape
=
3
*
self
.
nx
,
3
*
self
.
nx
self
.
dtype
=
np
.
complex
self
.
a_n
=
sample
.
dA
# other stuff for dense matrix computation
self
.
xyz
=
sample
.
xyz
self
.
belm
=
sample
.
belm
# self.belm = self.belm.reshape(len(self.belm) // 2, 2).T
# self.belm = self.belm.reshape(len(self.belm) // 2, 2).T
self
.
nv
=
sample
.
nv
self
.
pang
=
sample
.
pang
...
...
@@ -326,27 +338,23 @@ class DipolarOperator(LinearOperator):
# should read in/calculate dense matrix and not take it, will be done later
self
.
k
=
k
self
.
ksquared_mat
=
k
**
2
*
csr_matrix
(
diags
(
self
.
a_n
))
self
.
ksquared_mat
=
k
**
2
*
csr_matrix
(
diags
(
self
.
a_n
))
diag
=
self
.
a_n
.
copy
()
diag
[
self
.
boundary_nodes
]
=
np
.
zeros
(
self
.
nb
)
self
.
ksquared_matb
=
k
**
2
*
csr_matrix
(
diags
(
diag
))
self
.
ksquared_matb
=
k
**
2
*
csr_matrix
(
diags
(
diag
))
temp_dense
=
computedensematrix_kdep_py
(
self
.
belm
,
self
.
boundary_nodes
,
self
.
xyz
,
self
.
nv
,
self
.
pang
,
self
.
bndint_order
,
self
.
k
)
self
.
bndint_order
,
self
.
k
)
self
.
dense
=
np
.
reshape
(
temp_dense
,
(
self
.
nb
,
self
.
nb
))
def
_matvec
(
self
,
vec
):
if
self
.
k
==
None
:
raise
ValueError
(
"
No wave vector (k) has been assigned to the dynamic matrix yet. Do this by calling the set_k() method.
"
)
raise
ValueError
(
"
No wave vector (k) has been assigned to the dynamic matrix yet. Do this by calling the set_k() method.
"
)
else
:
return
N_dip
(
vec
,
self
.
k
,
self
.
ksquared_mat
,
self
.
ksquared_matb
,
self
.
nx
,
self
.
div_x
,
self
.
div_y
,
self
.
poisson
,
self
.
boundary_nodes
,
self
.
dense
,
self
.
laplace
,
self
.
grad_x
,
self
.
grad_y
,
self
.
a_n
)
self
.
div_x
,
self
.
div_y
,
self
.
poisson
,
self
.
boundary_nodes
,
self
.
dense
,
self
.
laplace
,
self
.
grad_x
,
self
.
grad_y
,
self
.
a_n
)
class
TotalDynMat
(
LinearOperator
):
...
...
@@ -366,7 +374,7 @@ class TotalDynMat(LinearOperator):
self
.
ksquared_mat
=
N_dip
.
ksquared_mat
self
.
nx
=
N_dip
.
nx
self
.
nb
=
N_dip
.
nb
self
.
shape
=
2
*
self
.
nx
,
2
*
self
.
nx
self
.
shape
=
2
*
self
.
nx
,
2
*
self
.
nx
self
.
dtype
=
np
.
complex
self
.
sparse_mat
=
N_nodip_mat
self
.
leftmulmat
=
leftmulmat
...
...
@@ -381,35 +389,36 @@ class TotalDynMat(LinearOperator):
self
.
nv
=
N_dip
.
nv
self
.
pang
=
N_dip
.
pang
super
().
__init__
(
self
.
dtype
,
self
.
shape
)
def
set_k
(
self
,
k
,
sparse_mat_k
):
# should read in/calculate dense matrix and not take it, will be done later
self
.
k
=
k
self
.
ksquared_mat
=
k
**
2
*
csr_matrix
(
diags
(
self
.
a_n
))
self
.
ksquared_mat
=
k
**
2
*
csr_matrix
(
diags
(
self
.
a_n
))
diag
=
self
.
a_n
.
copy
()
diag
[
self
.
boundary_nodes
]
=
np
.
zeros
(
self
.
nb
)
self
.
ksquared_matb
=
k
**
2
*
csr_matrix
(
diags
(
diag
))
self
.
ksquared_matb
=
k
**
2
*
csr_matrix
(
diags
(
diag
))
#self.ksquared_mat = k**2 * csr_matrix(identity(self.nx))
#self.dense = dense_k
#
self.ksquared_mat = k**2 * csr_matrix(identity(self.nx))
#
self.dense = dense_k
temp_dense
=
computedensematrix_kdep_py
(
self
.
belm
,
self
.
boundary_nodes
,
self
.
xyz
,
self
.
nv
,
self
.
pang
,
self
.
bndint_order
,
self
.
k
)
self
.
bndint_order
,
self
.
k
)
self
.
dense
=
np
.
reshape
(
temp_dense
,
(
self
.
nb
,
self
.
nb
))
self
.
sparse_mat_k
=
sparse_mat_k
def
_matvec
(
self
,
vec_loc
):
if
self
.
k
==
None
:
raise
ValueError
(
"
No wave vector (k) has been assigned to the dynamic matrix yet. Do this by calling the set_k() method.
"
)
raise
ValueError
(
"
No wave vector (k) has been assigned to the dynamic matrix yet. Do this by calling the set_k() method.
"
)
else
:
vec_lab
=
self
.
rightmulmat
.
dot
(
vec_loc
)
h_lab
=
self
.
sparse_mat_k
.
dot
(
vec_lab
)
+
N_dip
(
vec_lab
,
self
.
k
,
self
.
ksquared_mat
,
self
.
ksquared_matb
,
self
.
nx
,
self
.
div_x
,
self
.
div_y
,
self
.
poisson
,
self
.
boundary_nodes
,
self
.
dense
,
self
.
laplace
,
self
.
grad_x
,
self
.
grad_y
,
self
.
a_n
)
self
.
nx
,
self
.
div_x
,
self
.
div_y
,
self
.
poisson
,
self
.
boundary_nodes
,
self
.
dense
,
self
.
laplace
,
self
.
grad_x
,
self
.
grad_y
,
self
.
a_n
)
return
self
.
leftmulmat
.
dot
(
h_lab
)
...
...
@@ -417,6 +426,7 @@ class Error(Exception):
"""
Base class for exceptions in this module.
"""
pass
class
InvertError
(
Error
):
"""
Exception raised for errors in inversion of total dynamic matrix.
"""
...
...
@@ -432,7 +442,6 @@ class InvTotalDynMat(LinearOperator):
"""
Inverse of Dynamic Matrix
"""
def
__init__
(
self
,
demag_op
,
k
,
tol
=
1e-4
,
preconditioner
=
None
,
maxiter
=
1000
):
self
.
demag_op
=
demag_op
self
.
shape
=
demag_op
.
shape
self
.
dtype
=
np
.
complex
...
...
@@ -447,27 +456,27 @@ class InvTotalDynMat(LinearOperator):
logging
.
debug
(
"
inverse of {} requested
"
.
format
(
x
))
# reset demag_op so that it recalculates the demag field on first use
self
.
demag_op
.
counter
=
0
b
,
info
=
lgmres
(
self
.
demag_op
,
x
,
tol
=
self
.
tol
,
atol
=
1000
*
np
.
finfo
(
self
.
dtype
).
eps
,
b
,
info
=
lgmres
(
self
.
demag_op
,
x
,
tol
=
self
.
tol
,
atol
=
1000
*
np
.
finfo
(
self
.
dtype
).
eps
,
M
=
self
.
preconditioner
,
maxiter
=
self
.
maxiter
)
if
info
!=
0
:
raise
InvertError
(
self
.
k
,
info
,
self
.
counter
,
"
Error inverting dynamic matrix at k = {} : lgmres did not converge (info = {}).
"
"
The inverse was already calculated successfully {} times.
"
.
format
(
self
.
k
,
info
,
self
.
counter
))
raise
InvertError
(
self
.
k
,
info
,
self
.
counter
,
"
Error inverting dynamic matrix at k = {} : lgmres did not converge (info = {}).
"
"
The inverse was already calculated successfully {} times.
"
.
format
(
self
.
k
,
info
,
self
.
counter
))
self
.
counter
+=
1
#logging.info("inverted: {} times | gmres iterations: {}".format(self.counter, self.demag_op.counter))
#print("inverse computed")
#
logging.info("inverted: {} times | gmres iterations: {}".format(self.counter, self.demag_op.counter))
#
print("inverse computed")
return
b
class
SuperLUInv
(
LinearOperator
):
"""
SuperLU object as LinearOperator
"""
def
__init__
(
self
,
superlu
):
self
.
shape
=
superlu
.
shape
self
.
dtype
=
np
.
complex
self
.
superlu
=
superlu
def
_matvec
(
self
,
x
):
return
self
.
superlu
.
solve
(
x
)
...
...
@@ -476,14 +485,13 @@ def cross_operator(nx):
"""
Return the rotated cross product operator as a crs_matrix.
"""
col_idcs
=
np
.
empty
(
2
*
nx
,
dtype
=
np
.
int
)
row_idcs
=
np
.
empty
(
2
*
nx
,
dtype
=
np
.
int
)
data
=
np
.
empty
(
2
*
nx
)
col_idcs
=
np
.
empty
(
2
*
nx
,
dtype
=
np
.
int
)
row_idcs
=
np
.
empty
(
2
*
nx
,
dtype
=
np
.
int
)
data
=
np
.
empty
(
2
*
nx
)
data
[:
nx
]
=
-
1
row_idcs
[:
nx
]
=
np
.
arange
(
nx
)
col_idcs
[:
nx
]
=
nx
+
np
.
arange
(
nx
)
data
[
nx
:]
=
+
1
row_idcs
[
nx
:]
=
nx
+
np
.
arange
(
nx
)
col_idcs
[
nx
:]
=
np
.
arange
(
nx
)
return
csr_matrix
((
data
,
(
row_idcs
,
col_idcs
)),
shape
=
(
3
*
nx
,
3
*
nx
))
return
csr_matrix
((
data
,
(
row_idcs
,
col_idcs
)),
shape
=
(
3
*
nx
,
3
*
nx
))
This diff is collapsed.
Click to expand it.
src/tetrax/core/sample.py
+
3
−
0
View file @
e7769169
...
...
@@ -10,6 +10,7 @@ from ..helpers.math import normalize_single_3d_vec
class
Sample
():
def
__init__
(
self
,
name
=
"
my_sample
"
):
# material parameters
self
.
Msat
=
796e3
self
.
Aex
=
13e-12
self
.
name
=
name
...
...
@@ -22,12 +23,14 @@ class Sample():
self
.
Ku2
=
0.0
self
.
k_phi
=
0.0
self
.
k_theta
=
0.0
self
.
e_u
=
np
.
array
([
0
,
0
,
1
])
# uniaxial anistropy direction(s)
self
.
Kc1
=
0.0
self
.
v1Kc
=
np
.
array
([
1
,
0
,
0
])
self
.
v2Kc
=
np
.
array
([
0
,
1
,
0
])
self
.
v3Kc
=
np
.
array
([
0
,
0
,
1
])
# geometric parameters
self
.
geom
=
None
self
.
mesh
=
None
self
.
nx
=
None
...
...
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Save comment
Cancel
Please
register
or
sign in
to comment