Skip to content
Toggle navigation
P
Projects
G
Groups
S
Snippets
Help
Jigyasa Watwani
/
growth-pattern-control
This project
Loading...
Sign in
Toggle navigation
Go to a project
Project
Repository
Issues
0
Merge Requests
0
Pipelines
Wiki
Snippets
Members
Activity
Graph
Charts
Create a new issue
Jobs
Commits
Issue Boards
Files
Commits
Branches
Tags
Contributors
Graph
Compare
Charts
Commit
a2801d4f
authored
Aug 27, 2022
by
Jigyasa Watwani
Browse files
Options
Browse Files
Download
Email Patches
Plain Diff
numerical and exact solution
parent
4bb764e5
Show whitespace changes
Inline
Side-by-side
Showing
1 changed file
with
103 additions
and
0 deletions
moving_domain/moving_heat_equation_analytical.py
moving_domain/moving_heat_equation_analytical.py
0 → 100644
View file @
a2801d4f
import
numpy
as
np
import
matplotlib.pyplot
as
plt
from
matplotlib.widgets
import
Slider
import
progressbar
import
dolfin
as
df
df
.
set_log_level
(
df
.
LogLevel
.
ERROR
)
df
.
parameters
[
'form_compiler'
][
'optimize'
]
=
True
# parameters
k
=
1.0
T
=
1
dt
=
0.001
L
=
1
D
=
1.0
Nx
=
2000
Nt
=
1000
t
=
np
.
linspace
(
0
,
T
,
Nt
)
# diffusion and advection
def
diffusion
(
c
,
tc
):
return
(
D
*
df
.
inner
(
c
.
dx
(
0
),
tc
.
dx
(
0
)))
def
advection
(
c
,
tc
,
v
):
u
=
df
.
interpolate
(
v
,
c
.
function_space
())
return
(
df
.
inner
((
u
*
c
)
.
dx
(
0
),
tc
))
# create mesh
mesh
=
df
.
IntervalMesh
(
Nx
,
0
,
L
)
x
=
mesh
.
coordinates
()
# v = df.Constant(1.0)
v
=
df
.
Expression
(
'k*x[0]'
,
k
=
k
,
degree
=
1
)
# create function space
conc_element
=
df
.
FiniteElement
(
'P'
,
mesh
.
ufl_cell
(),
1
)
function_space
=
df
.
FunctionSpace
(
mesh
,
conc_element
)
# initial condition
c0
=
df
.
interpolate
(
df
.
Expression
(
'1 + 0.2*cos(pi*x[0]/L)'
,
pi
=
np
.
pi
,
L
=
L
,
degree
=
1
),
function_space
)
c0_array
=
c0
.
compute_vertex_values
(
mesh
)
# define variational problem
c
=
df
.
Function
(
function_space
)
tc
=
df
.
TestFunction
(
function_space
)
form
=
(
df
.
inner
((
c
-
c0
)
/
dt
,
tc
)
+
diffusion
(
c
,
tc
)
+
advection
(
c
,
tc
,
v
)
)
form
=
form
*
df
.
dx
# time stepping
ctot
=
np
.
zeros_like
(
t
)
x_array
=
np
.
zeros
((
len
(
t
),
mesh
.
num_vertices
()))
x_array
[
0
]
=
mesh
.
coordinates
()[:,
0
]
c_array
=
np
.
zeros
((
len
(
t
),
len
(
c0_array
)))
c_array
[
0
]
=
c0_array
ctot
[
0
]
=
df
.
assemble
(
c0
*
df
.
dx
(
mesh
))
for
n
in
progressbar
.
progressbar
(
range
(
1
,
len
(
t
))):
df
.
solve
(
form
==
0
,
c
)
c_array
[
n
]
=
c
.
compute_vertex_values
(
mesh
)
c0
.
assign
(
c
)
ctot
[
n
]
=
df
.
assemble
(
c0
*
df
.
dx
(
mesh
))
df
.
ALE
.
move
(
mesh
,
df
.
Expression
(
'v*dt'
,
v
=
v
,
dt
=
dt
,
degree
=
1
))
x_array
[
n
]
=
mesh
.
coordinates
()[:,
0
]
# analytical solution
c_exact
=
np
.
zeros
((
len
(
t
),
len
(
x
)))
for
i
in
range
(
len
(
t
)):
xprime
=
x_array
[
0
]
*
np
.
exp
(
-
k
*
t
[
i
])
/
L
tprime
=
(
D
/
(
2
*
k
*
L
**
2
))
*
(
1
-
np
.
exp
(
-
2
*
k
*
t
[
i
]))
# int = ((2*k**2*L**2)/(D**2))*(np.exp(-4*k*t[i])-1)
c_exact
[
i
]
=
np
.
exp
(
-
k
*
t
[
i
])
*
(
1
+
0.2
*
np
.
cos
(
np
.
pi
*
xprime
)
*
np
.
exp
(
-
np
.
pi
**
2
*
tprime
))
# plot c(x,t) computed numerically
fig
,
ax_comp
=
plt
.
subplots
(
1
,
1
,
figsize
=
(
8
,
6
))
ax_comp
.
set_xlabel
(
r'$x$'
)
ax_comp
.
set_ylabel
(
r'$c(x,t)$'
)
ax_comp
.
set_xlim
(
np
.
min
(
x_array
)
-
1
,
np
.
max
(
x_array
)
+
1
)
ax_comp
.
set_ylim
(
np
.
min
(
c_array
)
-
1
,
np
.
max
(
c_array
)
+
1
)
cplot
,
=
ax_comp
.
plot
(
x_array
[
0
],
c0_array
)
c_exactplot
,
=
ax_comp
.
plot
(
x_array
[
0
],
c_exact
[
0
],
'ro'
,
markevery
=
50
)
def
update
(
value
):
ti
=
np
.
abs
(
t
-
value
)
.
argmin
()
cplot
.
set_xdata
(
x_array
[
ti
])
cplot
.
set_ydata
(
c_array
[
ti
])
c_exactplot
.
set_xdata
(
x_array
[
ti
])
c_exactplot
.
set_ydata
(
c_exact
[
ti
])
plt
.
draw
()
sax
=
plt
.
axes
([
0.1
,
0.92
,
0.7
,
0.02
])
slider
=
Slider
(
sax
,
r'$t/\tau$'
,
min
(
t
),
max
(
t
),
valinit
=
min
(
t
),
valfmt
=
'
%3.1
f'
,
fc
=
'#999999'
)
slider
.
drawon
=
False
slider
.
on_changed
(
update
)
plt
.
show
()
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment