Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F86893090
helmholtz_problem_engine.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Oct 9, 06:14
Size
10 KB
Mime Type
text/x-python
Expires
Fri, Oct 11, 06:14 (2 d)
Engine
blob
Format
Raw Data
Handle
21501375
Attached To
R6746 RationalROMPy
helmholtz_problem_engine.py
View Options
# Copyright (C) 2018-2020 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
import
fenics
as
fen
from
.laplace_base_problem_engine
import
LaplaceBaseProblemEngine
from
rrompy.solver.fenics
import
fenZERO
,
fenONE
,
fenics2Sparse
from
rrompy.utilities.base
import
verbosityManager
as
vbMng
from
rrompy.utilities.exception_manager
import
RROMPyWarning
from
rrompy.parameter
import
parameterMap
as
pMap
__all__
=
[
'HelmholtzProblemEngine'
,
'ScatteringProblemEngine'
]
class
HelmholtzProblemEngine
(
LaplaceBaseProblemEngine
):
"""
Solver for generic Helmholtz problems with parametric wavenumber.
- \nabla \cdot (a \nabla u) - omega^2 * n**2 * u = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
cs: Numpy array representation of cs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
energyNormDualMatrix: Scipy sparse matrix representing dual inner
product between Riesz representers V-V.
degree_threshold: Threshold for ufl expression interpolation degree.
omega: Value of omega.
diffusivity: Value of a.
refractionIndex: Value of n.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
RobinDatumH: Value of h.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
dsToBeSet: Whether ds needs to be set.
"""
def
__init__
(
self
,
*
args
,
**
kwargs
):
super
()
.
__init__
(
*
args
,
**
kwargs
)
self
.
_affinePoly
=
True
self
.
nAs
=
2
self
.
parameterMap
=
pMap
([
2.
]
+
[
1.
]
*
(
self
.
npar
-
1
))
self
.
refractionIndex
=
fenONE
@property
def
refractionIndex
(
self
):
"""Value of n."""
return
self
.
_refractionIndex
@refractionIndex.setter
def
refractionIndex
(
self
,
refractionIndex
):
self
.
resetAs
()
if
not
isinstance
(
refractionIndex
,
(
list
,
tuple
,)):
refractionIndex
=
[
refractionIndex
,
fenZERO
]
self
.
_refractionIndex
=
refractionIndex
def
buildA
(
self
):
"""Build terms of operator of linear system."""
if
self
.
thAs
[
0
]
is
None
:
self
.
thAs
=
self
.
getMonomialWeights
(
self
.
nAs
)
if
self
.
As
[
0
]
is
None
:
self
.
autoSetDS
()
vbMng
(
self
,
"INIT"
,
"Assembling operator term A0."
,
20
)
DirichletBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZERO
,
self
.
DirichletBoundary
)
aRe
,
aIm
=
self
.
diffusivity
hRe
,
hIm
=
self
.
RobinDatumH
termNames
=
[
"diffusivity"
,
"RobinDatumH"
]
parsRe
=
self
.
iterReduceQuadratureDegree
(
zip
(
[
aRe
,
hRe
],
[
x
+
"Real"
for
x
in
termNames
]))
parsIm
=
self
.
iterReduceQuadratureDegree
(
zip
(
[
aIm
,
hIm
],
[
x
+
"Imag"
for
x
in
termNames
]))
a0Re
=
(
aRe
*
fen
.
dot
(
fen
.
grad
(
self
.
u
),
fen
.
grad
(
self
.
v
))
*
fen
.
dx
+
hRe
*
fen
.
dot
(
self
.
u
,
self
.
v
)
*
self
.
ds
(
1
))
a0Im
=
(
aIm
*
fen
.
dot
(
fen
.
grad
(
self
.
u
),
fen
.
grad
(
self
.
v
))
*
fen
.
dx
+
hIm
*
fen
.
dot
(
self
.
u
,
self
.
v
)
*
self
.
ds
(
1
))
self
.
As
[
0
]
=
(
fenics2Sparse
(
a0Re
,
parsRe
,
DirichletBC0
,
1
)
+
1.j
*
fenics2Sparse
(
a0Im
,
parsIm
,
DirichletBC0
,
0
))
vbMng
(
self
,
"DEL"
,
"Done assembling operator term."
,
20
)
if
self
.
As
[
1
]
is
None
:
vbMng
(
self
,
"INIT"
,
"Assembling operator term A1."
,
20
)
DirichletBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZERO
,
self
.
DirichletBoundary
)
nRe
,
nIm
=
self
.
refractionIndex
n2Re
,
n2Im
=
nRe
*
nRe
-
nIm
*
nIm
,
2
*
nRe
*
nIm
parsRe
=
self
.
iterReduceQuadratureDegree
(
zip
([
n2Re
],
[
"refractionIndexSquaredReal"
]))
parsIm
=
self
.
iterReduceQuadratureDegree
(
zip
([
n2Im
],
[
"refractionIndexSquaredImag"
]))
a1Re
=
-
n2Re
*
fen
.
dot
(
self
.
u
,
self
.
v
)
*
fen
.
dx
a1Im
=
-
n2Im
*
fen
.
dot
(
self
.
u
,
self
.
v
)
*
fen
.
dx
self
.
As
[
1
]
=
(
fenics2Sparse
(
a1Re
,
parsRe
,
DirichletBC0
,
0
)
+
1.j
*
fenics2Sparse
(
a1Im
,
parsIm
,
DirichletBC0
,
0
))
vbMng
(
self
,
"DEL"
,
"Done assembling operator term."
,
20
)
class
ScatteringProblemEngine
(
HelmholtzProblemEngine
):
"""
Solver for scattering problems with parametric wavenumber.
- \nabla \cdot (a \nabla u) - omega^2 * n**2 * u = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu +- i omega u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
cs: Numpy array representation of cs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
energyNormDualMatrix: Scipy sparse matrix representing dual inner
product between Riesz representers V-V.
degree_threshold: Threshold for ufl expression interpolation degree.
signR: Sign in ABC.
omega: Value of omega.
diffusivity: Value of a.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
dsToBeSet: Whether ds needs to be set.
"""
def
__init__
(
self
,
*
args
,
**
kwargs
):
self
.
silenceWarnings
=
True
super
()
.
__init__
(
*
args
,
**
kwargs
)
self
.
_affinePoly
=
True
del
self
.
silenceWarnings
self
.
nAs
=
3
self
.
parameterMap
=
pMap
(
1.
,
self
.
npar
)
self
.
signR
=
-
1.
@property
def
RobinDatumH
(
self
):
"""Value of h."""
if
not
hasattr
(
self
,
"silenceWarnings"
):
RROMPyWarning
(
"Scattering problems do not allow changes of h."
)
return
self
.
signR
@RobinDatumH.setter
def
RobinDatumH
(
self
,
RobinDatumH
):
if
not
hasattr
(
self
,
"silenceWarnings"
):
RROMPyWarning
((
"Scattering problems do not allow changes of h. "
"Ignoring assignment."
))
return
@property
def
signR
(
self
):
"""Value of signR."""
return
self
.
_signR
@signR.setter
def
signR
(
self
,
signR
):
self
.
resetAs
()
self
.
_signR
=
signR
def
buildA
(
self
):
"""Build terms of operator of linear system."""
if
self
.
thAs
[
0
]
is
None
:
self
.
thAs
=
self
.
getMonomialWeights
(
self
.
nAs
)
if
self
.
As
[
0
]
is
None
:
vbMng
(
self
,
"INIT"
,
"Assembling operator term A0."
,
20
)
DirichletBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZERO
,
self
.
DirichletBoundary
)
aRe
,
aIm
=
self
.
diffusivity
parsRe
=
self
.
iterReduceQuadratureDegree
(
zip
([
aRe
],
[
"diffusivityReal"
]))
parsIm
=
self
.
iterReduceQuadratureDegree
(
zip
([
aIm
],
[
"diffusivityImag"
]))
a0Re
=
aRe
*
fen
.
dot
(
fen
.
grad
(
self
.
u
),
fen
.
grad
(
self
.
v
))
*
fen
.
dx
a0Im
=
aIm
*
fen
.
dot
(
fen
.
grad
(
self
.
u
),
fen
.
grad
(
self
.
v
))
*
fen
.
dx
self
.
As
[
0
]
=
(
fenics2Sparse
(
a0Re
,
parsRe
,
DirichletBC0
,
1
)
+
1.j
*
fenics2Sparse
(
a0Im
,
parsIm
,
DirichletBC0
,
0
))
vbMng
(
self
,
"DEL"
,
"Done assembling operator term."
,
20
)
if
self
.
As
[
1
]
is
None
:
self
.
autoSetDS
()
vbMng
(
self
,
"INIT"
,
"Assembling operator term A1."
,
20
)
DirichletBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZERO
,
self
.
DirichletBoundary
)
a1
=
fen
.
dot
(
self
.
u
,
self
.
v
)
*
self
.
ds
(
1
)
self
.
As
[
1
]
=
(
self
.
signR
*
1.j
*
fenics2Sparse
(
a1
,
{},
DirichletBC0
,
0
))
vbMng
(
self
,
"DEL"
,
"Done assembling operator term."
,
20
)
if
self
.
As
[
2
]
is
None
:
vbMng
(
self
,
"INIT"
,
"Assembling operator term A2."
,
20
)
DirichletBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZERO
,
self
.
DirichletBoundary
)
nRe
,
nIm
=
self
.
refractionIndex
n2Re
,
n2Im
=
nRe
*
nRe
-
nIm
*
nIm
,
2
*
nRe
*
nIm
parsRe
=
self
.
iterReduceQuadratureDegree
(
zip
([
n2Re
],
[
"refractionIndexSquaredReal"
]))
parsIm
=
self
.
iterReduceQuadratureDegree
(
zip
([
n2Im
],
[
"refractionIndexSquaredImag"
]))
a2Re
=
-
n2Re
*
fen
.
dot
(
self
.
u
,
self
.
v
)
*
fen
.
dx
a2Im
=
-
n2Im
*
fen
.
dot
(
self
.
u
,
self
.
v
)
*
fen
.
dx
self
.
As
[
2
]
=
(
fenics2Sparse
(
a2Re
,
parsRe
,
DirichletBC0
,
0
)
+
1.j
*
fenics2Sparse
(
a2Im
,
parsIm
,
DirichletBC0
,
0
))
vbMng
(
self
,
"DEL"
,
"Done assembling operator term."
,
20
)
Event Timeline
Log In to Comment