Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 44db013c authored by SAMASSA Karamoko's avatar SAMASSA Karamoko
Browse files

fixing sympy expression creation from a string

parent 9fda87d3
No related branches found
No related tags found
1 merge request!10Draft: Feature/upgrade to fenicsx
Hide whitespace changes
Inline Side-by-side
src/bvpy/domains/abstract.py
+ 1
− 0
View file @ 44db013c
......@@ -403,6 +403,7 @@ class AbstractDomain(ABC):
self.model.addPhysicalGroup(dim, [self.volumes[0]])
self.model.mesh.generate(dim)
self.mesh, cell_markers, facet_markers = gmshio.model_to_mesh(self.model, MPI.COMM_SELF, 0)
# self.bdata = _generate_default_boundary_data(self.mesh)
self.ds = ufl.Measure('ds', domain=self.mesh)
......
src/bvpy/ibvp.py
+ 20
− 20
View file @ 44db013c
......@@ -75,7 +75,7 @@ class IBVP(BVP):
self.time = Time()
self.scheme = scheme(vform)
self.previous_solution = self.solution.copy(deepcopy=True)
self.previous_solution = self.solution.copy()
self.scheme.set_previous_solution(self.previous_solution)
self._progressbar = ProgressBar()
self.solution_steps = None
......@@ -112,31 +112,31 @@ class IBVP(BVP):
if store_steps:
self.solution_steps = OrderedDict()
self.solution_steps[self.time.current] =\
self.solution.copy(deepcopy=True)
self.solution.copy()
if filename is not None:
# file = fe.XDMFFile(MPI_COMM_WORLD, filename)
file = io.XDMFFile(MPI_COMM_WORLD, filename)
file.parameters['functions_share_mesh'] = True
file.write(self.solution, self.time.current)
logger.setLevel(logging.WARNING)
self._progressbar.set_range(t_i, t_f)
while t_f-self.time.current > dt:
self.step(dt, **kwargs)
file.write(self.solution, self.time.current)
if store_steps:
self.solution_steps[self.time.current] =\
self.solution.copy(deepcopy=True)
with io.XDMFFile(MPI_COMM_WORLD, filename, "w") as file:
file.write_mesh(self.solution.function_space.mesh)
# file.parameters['functions_share_mesh'] = True
file.write_function(self.solution, self.time.current)
logger.setLevel(logging.WARNING)
self._progressbar.set_range(t_i, t_f)
while t_f-self.time.current > dt:
self.step(dt, **kwargs)
file.write_function(self.solution, self.time.current)
if store_steps:
self.solution_steps[self.time.current] =\
self.solution.copy()
self._progressbar.update(self.time.current)
self._progressbar.show()
self.step(t_f-self.time.current, **kwargs)
file.write_function(self.solution, self.time.current)
self._progressbar.update(self.time.current)
self._progressbar.show()
self.step(t_f-self.time.current, **kwargs)
file.write(self.solution, self.time.current)
self._progressbar.update(self.time.current)
self._progressbar.show()
file.close()
else:
logger.setLevel(logging.WARNING)
self._progressbar.set_range(t_i, t_f)
......
src/bvpy/solvers/linear.py
+ 1
− 1
View file @ 44db013c
......@@ -96,7 +96,7 @@ class LinearSolver():
logger.info(' Preconditioner: '
+ str(self.parameters['preconditioner']))
logger.info(' Size : '
+ str(self._u.function_space().dim()))
+ str(self._u.function_space.mesh.geometry.dim))
logger.info(' Solving LinearProblem ...')
# A = fe.assemble(self._a)
......
src/bvpy/solvers/time_integration.py
+ 12
− 8
View file @ 44db013c
......@@ -22,6 +22,7 @@
# import fenics as fe
from ufl.classes import Zero
import ufl
from dolfinx import fem
class Time():
"""Container class for Time.
......@@ -104,19 +105,22 @@ class ImplicitEuler(FirstOrderScheme):
def apply(self, u, v, sol):
lhs, rhs = self._vform.get_form(u, v, sol)
if isinstance(rhs, Zero):
dudt = ufl.dot((sol-self.previous_solution)
/ ufl.Constant(self.dt), v) * self.dx
# dudt = ufl.dot((sol-self.previous_solution)
# / ufl.Constant(self.dt), v) * self.dx
dudt = ufl.dot((sol-self.previous_solution)/self.dt, v) * self.dx
lhs += dudt
else:
dudt = ufl.dot((u-self.previous_solution)
/ ufl.Constant(self.dt), v) * self.dx
# dudt = ufl.dot((u-self.previous_solution)
# / ufl.Constant(self.dt), v) * self.dx
dudt = ufl.dot((u-self.previous_solution)/self.dt, v) * self.dx
dudt_lhs, dudt_rhs = ufl.lhs(dudt), ufl.rhs(dudt)
lhs += dudt_lhs
rhs += dudt_rhs
dudt_lhs, dudt_rhs = fem.form(ufl.lhs(dudt)), fem.form(ufl.rhs(dudt))
# lhs += dudt_lhs
# rhs += dudt_rhs
lhs = dudt_lhs
rhs = dudt_rhs
return lhs, rhs
......
src/bvpy/utils/pre_processing.py
+ 5
− 3
View file @ 44db013c
......@@ -24,7 +24,8 @@ import logging
# import fenics as fe
import numpy as np
from sympy import Function, Symbol, lambdify, symbols, sympify
import sympy
from sympy import Function, Symbol
from sympy.parsing.sympy_parser import parse_expr
from types import FunctionType
......@@ -45,8 +46,9 @@ def _expression_from_string(source, functionspace, degree, **kwargs):
for src in source:
expr.append(str(parse_expr(src, local_dict=_XYZ, evaluate=False)))
else:
# expr = str(parse_expr(source, local_dict=_XYZ, evaluate=False))
expr = sympify(source)
expr = parse_expr(source, local_dict=_XYZ, evaluate=False)
x = Symbol('x')
expr = sympy.lambdify(x, expr)
#TODO: is this part useful anymore ?
......
tutorials/bvpy_tutorial_5_reaction_diffusion.ipynb
+ 4
− 3
View file @ 44db013c
%% Cell type:markdown id: tags:
# Reaction-diffusion systems
This tutorial tackles a classic and yet complex example of *bvp* highly influencial in mathematical biology: Reaction-diffusion systems. These encompases any *PDE* system featuring a Laplacian (the diffusion term) and a source (the reaction term), fonction of the unknown field(s). Their general form reads:
$$
\partial_t\phi = D\Delta\phi + f(\phi)
$$
The reaction term $f(\phi)$ can be highly non-linear. As a matter of fact, non-linear reaction terms can generate very complex dynamics and non-trivial stationnary solutions.
Through this tutorial we want to show how **bvpy** can be used to simulate such equations. In particular we are going to implement two classic reaction-diffusion system: A wave-front propagation and a Turing pattern.
**Covered topics:**
- The `TransportForm` an `CoupledTransportForm` classes from the `bvpy.vforms` module.
- Basic use of the `ibvp` class to implement initial-boundary-value problems.
- Basic use of the `SolutionExplorer` class from the `bvpy.utils.post_processing` sub-module to handle simulation results.
Download the notebook: [bvpy_tutorial_5](https://gitlab.inria.fr/mosaic/publications/bvpy/-/raw/develop/tutorials/bvpy_tutorial_5_turing_system.ipynb?inline=false).
%% Cell type:markdown id: tags:
## Wave-front propagation
Let's consider the following reaction-diffusion equation:
$$
\begin{array}{c}
\partial_t \phi = D\Delta\phi + f(\phi) \\
\text{with:} \quad\quad f(\phi) = -\displaystyle\prod_{i=0}^{2}(\phi-\phi_i)
\end{array}
$$
For a deep understanding of the importance of this type of equation in biology and biochemistry, the interested reader is directed towards the books **Mathematical Biology I & II** by *J.D. Murray*.
%% Cell type:markdown id: tags:
### Problem implementation
%% Cell type:markdown id: tags:
**Domain:** Let's consider a simple continuous flat band of aspect ration 5 in the x direction:
%% Cell type:code id: tags:
``` python
from bvpy.domains import Rectangle
domain = Rectangle(length=1, width=.2, cell_size=0.02, dim=2, clear=True)
```
%% Cell type:markdown id: tags:
**Vform:** To implement the reaction-diffusion equation, we will use the `TransportForm` class, where we will simply set a value for the diffusion coefficient and implement the polynomial expression of the reaction function:
%% Cell type:code id: tags:
``` python
from bvpy.vforms import TransportForm
reaction_diffusion = TransportForm(diffusion_coef=.01, reaction= lambda u: -u*(u-.2)*(u-1))
```
%% Cell type:markdown id: tags:
**Boundary conditions:**
%% Cell type:code id: tags:
``` python
from bvpy.boundary_conditions.dirichlet import dirichlet
fixed_boundary_values = [dirichlet(1, "near(x, 0)"),
dirichlet(0, "near(x, 1)")]
```
%% Cell type:markdown id: tags:
**Initial condition:** Since time is an integration variable in *ibvps*, we need to set an initial condition. Since we want to simulate the propagation of an interface between two regions
%% Cell type:code id: tags:
``` python
from bvpy.utils.pre_processing import create_expression
import sympy
init = create_expression("1+tanh(-x/.1)", degree=4)
```
%% Cell type:markdown id: tags:
**IBVP implementation:**
%% Cell type:code id: tags:
``` python
from bvpy import IBVP
wave_propagation = IBVP(domain, reaction_diffusion, bc=fixed_boundary_values, initial_state=init)
```
%% Cell type:markdown id: tags:
### Time integration
Once the problem is set, we can integrate it over a time intervalle $\Delta t = [t_{min}, t_{max}]$ with a time resolution $dt$:
%% Cell type:code id: tags:
``` python
t_min = 0
t_max = 29
dt = 0.1
wave_propagation.integrate(t_min, t_max, dt, '../data/tutorial_results/tuto_5_result_1.xdmf', store_steps=True)
```
%% Cell type:markdown id: tags:
>**Notes:**
> - For `ibvps` saving the solution with the `save()` function is not required anymore for recording is made *de facto* by the `.integrate()` method.
> - When the `store_steps` argument is set to `True`, all integration steps are recorded, otherwise only the final state is.
%% Cell type:markdown id: tags:
### Solution visualization
%% Cell type:markdown id: tags:
Visualization of temporal dynamics is not possible in the current version of **Bvpy**, we therefore need an external software to assess the qualitatively the result of the simulation. This can be done, for instance, with the help of the **Paraview** software.
%% Cell type:code id: tags:
``` python
from IPython.display import Image
Image(filename='../data/tutorial_results/tuto_5_result_1_paraview_snapshot.gif.png')
```
%% Cell type:markdown id: tags:
The animated gif above shows the time evolution of the solution $\phi\colon(t,\mathbf{x})\to\phi(t,\mathbf{x})$ of the `wave_propagation` problem above. The heatmap encodes the solution amplitude: red for 1 and blue for 0.
%% Cell type:markdown id: tags:
More quantitative visual representation can be achieved with tools from our library. For instance, one might be interesteed in ploting the amplitude profiles $\phi(t_i, \mathbf{x}(\xi))$ for various time steps $t_i$ and with $\xi$ the curvilinear abscisse along a specific curved within the considered domain. For instance, we can plot the $\phi$ profile at various times along the horizontal axe, *i.e.* the family of curves:
$$
\begin{array}{rl}
\phi_i\colon & [0, 1] \mapsto [0, 1] \\
& x \to \phi(t_i, [x, 0])
\end{array}
$$
%% Cell type:markdown id: tags:
Solution analysis can be done with the help of the `SolutionExplorer` class from the `bvpy.utils.post_processing` sub-module. This class roughly wraps the FEniCS technicity and enable a series of basic manipulations in a intuitive way.
One of its useful method is the `.geometric_filter()` one. from a simple geometric rule on the domain (*e.g.* "x>=0" ) it returns an array of the solution values on all the vertices contained within the corresponding sub-domain. For instance, here below, the geometric filter returns the solution on all vertices on the line "y=0". The `sort="x"` argument insures that the list is sorted by increasing x value.
%% Cell type:code id: tags:
``` python
from bvpy.utils.post_processing import SolutionExplorer
concentration_over_time = []
for concentration in wave_propagation.solution_steps.values():
concentration_over_time.append(SolutionExplorer(concentration).geometric_filter(axe="y",
threshold=0,
comparator="==",
sort="x"))
```
%% Cell type:markdown id: tags:
Once this projection over the abscissae axe is done for each time step, one can plot the resulting list of arrays:
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(15, 5))
for i, c in enumerate(concentration_over_time[::30]):
plt.plot(c, marker="o", label=f"{30*i}th")
plt.xlabel('Abscissae (A.U.)')
plt.ylabel('Concentration (A.U.)')
plt.title('Concentration horizontal distribution over time.', fontsize=16)
fig.legend(loc='right', title="time steps:")
fig.show()
```
%% Cell type:markdown id: tags:
## Turing patterns
Let's now try to implement another iconic example of *ibvp* in mathematical biology: Turing reaction-diffusion systems.
%% Cell type:markdown id: tags:
Such system are usually described as a coupled set of two *PDEs* depicting the spatio-temporal dynamics of two scalar fields. Their general form reads:
$$
\begin{cases}
\partial_t a = D_a \Delta a + r_a(a, b) \\
\partial_t b = D_b \Delta b + r_b(a, b) \\
\end{cases}
$$
where the *reaction terms* $r_a(a,b)$ and $r_b(a,b)$ are usually highly non-linear smooth functions.
In chemistry, biochemistry and developmental biology such systems can formalize chemical componds diffusing and reacting with each other, typically in a catalytic manner.
In the present example, we are considering one specific case of Turing system: The **Lengyel-Epstein model** [(Lengyel & Epstein, Science (1991))](https://dx.doi.org/10.1126/science.251.4994.650) which has been recently used to investigate the zebra pattern formation [(Jeong *et al*, Statistical Mechanics and its Applications (2017))](https://dx.doi.org/10.1016/j.physa.2017.02.014).
The particularity of the **Lengyel-Epstein model** is to use very specific forms for the reaction terms:
$$
r_a(a,b)=k_a b\Big(1-\frac{a}{1+b^2}\Big) \quad\quad \text{and} \quad\quad r_b(a,b) = k_b -b -\frac{4a}{1+b^2}.
\label{eq:LE_model_source_terms}
$$
%% Cell type:markdown id: tags:
### Domain definition and boundary conditions
Let's consider a simple square domain:
%% Cell type:code id: tags:
``` python
from bvpy.domains import Rectangle
skin_patch = Rectangle(length=10, width=10, cell_size=0.2, dim=2, clear=True)
```
%% Cell type:markdown id: tags:
But with periodic boundary conditions to ease the pattern formation.
Such boundary conditions can be implemented thanks to the `SquarePeriodicBoundary` from the `bvpy.boundary_conditions.periodic` sub-module:
%% Cell type:code id: tags:
``` python
from bvpy.boundary_conditions.periodic import SquarePeriodicBoundary
periodic_bc=SquarePeriodicBoundary(length=10, width=10)
```
%% Cell type:markdown id: tags:
### Variational formulation of the Lengyel-Epstein model
%% Cell type:markdown id: tags:
As demonstrated in the previous section, reaction-diffusion problems can be implemented with the `TransportForm`. But in the present case, we are interested in the system of two reaction-diffusion equation that are coupled to each other. We derived a specific class to deal with such problems: The `CoupledTransportForm` from the `bvpy.vforms.transport` sub-module:
%% Cell type:code id: tags:
``` python
from bvpy.vforms import CoupledTransportForm
LE_model = CoupledTransportForm()
LE_model.info()
```
%% Cell type:markdown id: tags:
From the `.info()` method above, we get that we need to define three things:
- A diffusion coefficient.
- A reaction term.
- A source term.
For the reaction terms we will use the functions given in expression ($\ref{eq:LE_model_source_terms}$) and we will take a null source term.
> **Note:** The difference that is usually made between a *source* and a *reaction* term is that the former is either constant or position dependent on the considered domain whereas the latter is explicitly a function of the unkown field, *i.e.* $s=s(\mathbf{x})$ and $r=r(a,b)$.
The particularity of the `CoupledTransportForm` is that it features a `.add_species()` method that enables to add as much equations (and the corresponding unknown fields) as desired. The arguments of this methods are:
- The name of the unknown.
- All the coefficients (*i.e.* diffusion coefficient, soure and reaction terms) of the corresponding equations that need to be set.
%% Cell type:code id: tags:
``` python
Da, Db = 1, .02
ka, kb = 9, 11
LE_model.add_species('a', Da, lambda a, b: ka*(b-(a*b)/(1+b*b)))
LE_model.add_species('b', Db, lambda a, b: kb-b-(4*a*b)/(1+b*b))
```
%% Cell type:markdown id: tags:
> **Notes:**
> - As shown above, we can pass any function that we want to define the required coefficients.
> - The source term is set to 0 by default, see method documentation.
%% Cell type:markdown id: tags:
### Initial conditions
We will start here from two random distributions $a_0(\mathbf{x})$ and $b_0(\mathbf{x})$ for the unknown fields.
%% Cell type:code id: tags:
``` python
import numpy as np
from bvpy.utils.pre_processing import create_expression
np.random.seed(1093)
skin_patch.discretize()
random = np.random.uniform(low=-1, high=1, size=skin_patch.mesh.num_entities(2))
random_pattern = [[1+0.04*ka**2+0.1*r, 0.2*kb+0.1*r] for r in random]
initial_pattern = create_expression(random_pattern)
```
%% Cell type:markdown id: tags:
> **Note:** to generate the random pattern we need a grid of points. These points correspond to the vertices of the meshed domain. But at this stage, the domain is not meshed yet. This meshing procedure is performed by the `.discretize()` method.
%% Cell type:markdown id: tags:
### Problem instanciation and integration
%% Cell type:code id: tags:
``` python
from bvpy.ibvp import IBVP
from bvpy.solvers.time_integration import ImplicitEuler
turing_system = IBVP(skin_patch, LE_model, periodic_bc, initial_pattern)
```
%% Cell type:code id: tags:
``` python
t_min = 0
t_max = 9
dt = 0.1
save_path = '../data/tutorial_results/tuto_5_result_2.xdmf'
turing_system.integrate(t_min, t_max, dt, save_path)
```
%% Cell type:markdown id: tags:
> **Note:** With *ibvps* no need to call the `save()` function from the `bvpy.utils.io` sub-module for the `.integrate()` method can save the results directly.
%% Cell type:markdown id: tags:
### Results visualization and analysis
Let's have a look at the final patterns.
%% Cell type:code id: tags:
``` python
from bvpy.utils.visu import plot
final_pattern_a = turing_system.solution.sub(0)
final_pattern_b = turing_system.solution.sub(1)
plot(final_pattern_a)
```
%% Cell type:code id: tags:
``` python
plot(final_pattern_b, colorscale="viridis")
```
%% Cell type:markdown id: tags:
Let's have a look at the whole dynamics, recorded as a `gif` within **Paraview**:
%% Cell type:code id: tags:
``` python
from IPython.display import Image
Image(filename='../data/tutorial_results/tuto_5_result_2_paraview_snapshot.gif.png')
#from IPython.display import HTML
#HTML('<img src="../data/tutorial_results/tuto_5_result_2_paraview_snapshot.gif">')
```
......
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