Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F61683710
laplace_disk_gaussian.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, May 8, 07:44
Size
7 KB
Mime Type
text/x-python
Expires
Fri, May 10, 07:44 (2 d)
Engine
blob
Format
Raw Data
Handle
17544627
Attached To
R6746 RationalROMPy
laplace_disk_gaussian.py
View Options
# Copyright (C) 2018 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
numpy
as
np
import
fenics
as
fen
from
rrompy.utilities.base.types
import
Np1D
,
Tuple
,
FenExpr
,
paramVal
from
.laplace_base_problem_engine
import
LaplaceBaseProblemEngine
from
rrompy.solver.fenics
import
fenZERO
,
fenONE
from
rrompy.utilities.base
import
verbosityDepth
from
rrompy.utilities.poly_fitting.polynomial
import
(
hashDerivativeToIdx
as
hashD
,
hashIdxToDerivative
as
hashI
)
from
rrompy.solver.fenics
import
fenics2Vector
__all__
=
[
'LaplaceDiskGaussian'
]
class
LaplaceDiskGaussian
(
LaplaceBaseProblemEngine
):
"""
Solver for disk Laplace problems with parametric forcing term center.
- \Delta u = C exp(-.5 * ||\cdot - (mu, 0)||^2) in \Omega = B(0, 5)
u = 0 on \partial\Omega.
"""
def
__init__
(
self
,
n
:
int
,
mu0
:
paramVal
=
[
0.
],
degree_threshold
:
int
=
np
.
inf
,
verbosity
:
int
=
10
,
timestamp
:
bool
=
True
):
super
()
.
__init__
(
mu0
=
mu0
,
degree_threshold
=
degree_threshold
,
verbosity
=
verbosity
,
timestamp
=
timestamp
)
self
.
nbs
=
20
self
.
computebsFactors
()
self
.
forcingTermMu
=
np
.
nan
import
mshr
mesh
=
mshr
.
generate_mesh
(
mshr
.
Circle
(
fen
.
Point
(
0.
,
0.
),
5.
),
3
*
n
)
self
.
V
=
fen
.
FunctionSpace
(
mesh
,
"P"
,
1
)
def
getForcingTerm
(
self
,
mu
:
paramVal
=
[])
->
Tuple
[
FenExpr
,
FenExpr
]:
"""Compute forcing term."""
mu
=
self
.
checkParameter
(
mu
)
if
mu
!=
self
.
forcingTermMu
:
if
self
.
verbosity
>=
25
:
verbosityDepth
(
"INIT"
,
(
"Assembling base expression for "
"forcing term."
),
timestamp
=
self
.
timestamp
)
x
,
y
=
fen
.
SpatialCoordinate
(
self
.
V
.
mesh
())[:]
C
=
np
.
exp
(
-.
5
*
mu
(
0
,
0
)
**
2.
)
CR
,
CI
=
np
.
real
(
C
),
np
.
imag
(
C
)
f0
=
(
2
*
np
.
pi
)
**
-.
5
*
fen
.
exp
(
-.
5
*
(
x
**
2.
+
y
**
2.
))
muR
,
muI
=
np
.
real
(
mu
(
0
,
0
)),
np
.
imag
(
mu
(
0
,
0
))
f1R
=
fen
.
exp
(
muR
*
x
)
*
fen
.
cos
(
muI
*
x
)
f1I
=
fen
.
exp
(
muR
*
x
)
*
fen
.
sin
(
muI
*
x
)
self
.
forcingTerm
=
[
f0
*
(
CR
*
f1R
-
CI
*
f1I
),
f0
*
(
CR
*
f1I
+
CI
*
f1R
)]
self
.
forcingTermMu
=
mu
if
self
.
verbosity
>=
25
:
verbosityDepth
(
"DEL"
,
"Done assembling base expression."
,
timestamp
=
self
.
timestamp
)
return
self
.
forcingTerm
def
computebsFactors
(
self
):
self
.
bsFactors
=
np
.
zeros
((
self
.
nbs
,
self
.
nbs
),
dtype
=
float
)
self
.
bsFactors
[
0
,
0
]
=
1.
self
.
bsFactors
[
1
,
1
]
=
1.
for
j
in
range
(
2
,
self
.
nbs
):
l
=
(
j
+
1
)
%
2
+
1
J
=
np
.
arange
(
l
,
j
+
1
,
2
)
self
.
bsFactors
[
j
,
J
]
=
self
.
bsFactors
[
j
-
1
,
J
-
1
]
if
l
==
2
:
l
=
0
J
=
np
.
arange
(
l
,
j
,
2
)
self
.
bsFactors
[
j
,
J
]
+=
np
.
multiply
(
-
1
-
J
,
self
.
bsFactors
[
j
-
1
,
J
+
1
])
self
.
bsFactors
[
j
,
l
:
j
+
2
:
2
]
/=
j
def
getExtraFactorB
(
self
,
mu
:
paramVal
=
[],
derI
:
int
=
0
)
->
Tuple
[
FenExpr
,
FenExpr
]:
"""Compute extra expression in RHS."""
mu
=
self
.
checkParameter
(
mu
)
if
self
.
verbosity
>=
25
:
verbosityDepth
(
"INIT"
,
(
"Assembling auxiliary expression for "
"forcing term derivative."
),
timestamp
=
self
.
timestamp
)
muR
,
muI
=
np
.
real
(
mu
(
0
,
0
)),
np
.
imag
(
mu
(
0
,
0
))
x
=
fen
.
SpatialCoordinate
(
self
.
V
.
mesh
())[
0
]
l
=
derI
%
2
if
l
==
0
:
powR
,
powI
=
fenONE
,
fenZERO
else
:
powR
,
powI
=
x
-
muR
,
fen
.
Constant
(
muI
)
exprR
,
exprI
=
[
self
.
bsFactors
[
derI
,
l
]
*
k
for
k
in
[
powR
,
powI
]]
for
j
in
range
(
l
+
2
,
derI
+
1
,
2
):
for
_
in
range
(
2
):
powR
,
powI
=
(
powR
*
(
x
-
muR
)
-
powI
*
muI
,
powR
*
muI
+
powI
*
(
x
-
muR
))
exprR
+=
self
.
bsFactors
[
derI
,
j
]
*
powR
exprI
+=
self
.
bsFactors
[
derI
,
j
]
*
powI
if
self
.
verbosity
>=
25
:
verbosityDepth
(
"DEL"
,
"Done assembling auxiliary expression."
,
timestamp
=
self
.
timestamp
)
return
[
exprR
,
exprI
]
def
b
(
self
,
mu
:
paramVal
=
[],
der
:
int
=
0
,
homogeneized
:
bool
=
False
)
->
Np1D
:
"""Assemble (derivative of) RHS of linear system."""
mu
=
self
.
checkParameter
(
mu
)
if
not
hasattr
(
der
,
"__len__"
):
der
=
[
der
]
*
self
.
npar
derI
=
hashD
(
der
)
nbsTot
=
self
.
nbsH
if
homogeneized
else
self
.
nbs
bs
=
self
.
bsH
if
homogeneized
else
self
.
bs
if
homogeneized
and
self
.
mu0
!=
self
.
mu0BC
:
self
.
liftDirichletData
(
self
.
mu0
)
for
j
in
range
(
derI
,
nbsTot
):
if
bs
[
j
]
is
None
:
if
self
.
verbosity
>=
20
:
verbosityDepth
(
"INIT"
,
(
"Assembling forcing term "
"b{}."
)
.
format
(
j
),
timestamp
=
self
.
timestamp
)
if
j
<
self
.
nbs
:
fRe
,
fIm
=
self
.
getForcingTerm
(
self
.
mu0
)
cRe
,
cIm
=
self
.
getExtraFactorB
(
self
.
mu0
,
j
)
cfRe
,
cfIm
=
cRe
*
fRe
-
cIm
*
fIm
,
cRe
*
fIm
+
cIm
*
fRe
else
:
cfRe
,
cfIm
=
fenZERO
,
fenZERO
parsRe
=
self
.
iterReduceQuadratureDegree
(
zip
([
cfRe
],
[
"forcingTermDer{}Real"
.
format
(
j
)]))
parsIm
=
self
.
iterReduceQuadratureDegree
(
zip
([
cfIm
],
[
"forcingTermDer{}Imag"
.
format
(
j
)]))
L0Re
=
fen
.
dot
(
cfRe
,
self
.
v
)
*
fen
.
dx
L0Im
=
fen
.
dot
(
cfIm
,
self
.
v
)
*
fen
.
dx
DirichletBC0
=
fen
.
DirichletBC
(
self
.
V
,
fenZERO
,
self
.
DirichletBoundary
)
b
=
(
fenics2Vector
(
L0Re
,
parsRe
,
DirichletBC0
,
1
)
+
1.j
*
fenics2Vector
(
L0Im
,
parsIm
,
DirichletBC0
,
1
))
if
homogeneized
:
Ader
=
self
.
A
(
self
.
mu0
,
hashI
(
j
,
self
.
npar
))
b
-=
Ader
.
dot
(
self
.
liftedDirichletDatum
)
if
homogeneized
:
self
.
bsH
[
j
]
=
b
else
:
self
.
bs
[
j
]
=
b
if
self
.
verbosity
>=
20
:
verbosityDepth
(
"DEL"
,
"Done assembling forcing term."
,
timestamp
=
self
.
timestamp
)
return
self
.
_assembleb
(
mu
,
der
,
derI
,
homogeneized
,
self
.
mu0
)
Event Timeline
Log In to Comment