Defined functions in initial conditions and equation's coefficient

[DanielaScarpelli]

Dear Dr. Zwicker,

I am really grateful on your develpment on py-pde. I was trying to use your package to solve the following PDE:

The coefficient \sigma depends on the coordinate r and is written in terms of many other defined (def) functions (temperature, density, etc..). The same goes for the velocity v. The spatial coordinate r goes from R_i to R (this provides the grid)

The initial condition is given by a piecewise function:

The boundary conditions are:

1. s(R_i,t) = 0

My questions are:

a) Is it possible with py-pde to specify an initial condition as the one I have, that is, given by a piecewise function? If so, how?

b) Does py-pde solve a PDE as mine where the coefficients (in these case \sigma and v) are defined in terms of other functions. If so, how could it be implemented?

I have searched, without success, for these answers on the several examples provided and in the user and reference manuals. I would very much appreciate your help.

Looking forward to your answer.

Dr. Daniela Pérez

PS: the pde above models the evolution of the magnetic field in neutron stars.

[david-zwicker]

Dear Daniela,

yes, both aspects should be possible!

a) I just updated the package (so please use the latest version from GitHub for now), so you can initialize a field from a piecewise expression. The syntax is determined by sympy:

grid = CartesianGrid([[0, 4]], 32)
field = ScalarField.from_expression(grid, "Piecewise((x**2, x>2), (1+x, x<=2))")

The field will now be linear in x for x<2 and quadratic in x for larger values.

b) Coefficients that depend on field values are certainly possible. You can simply use mathematical expressions in the definitions of general PDEs using the PDE class. If you need more flexibility, you might need to implement your own model class by subclassing PDEBase. The examples might be a good starting point for both things.

[DanielaScarpelli]

Dear David,

thanks a lor for your prompt reply! I will then follow the syntax in a) and try to implement your suggestion for b).

Best,

Daniela

[DanielaScarpelli]

Dear David,

I tried to implement your suggestion:

grid = CartesianGrid([[0, 4]], 32)
field = ScalarField.from_expression(grid, "Piecewise((x**2, x>2), (1+x, x<=2))")

But I get the following error:

Traceback (most recent call last):
File "/home/dperez/PycharmProjects/diffusion_equation/main.py", line 13, in
field = ScalarField.from_expression(grid, "Piecewise((x**2, x>2), (1+x, x<=2))")
File "/home/dperez/.local/lib/python3.8/site-packages/pde/fields/scalar.py", line 61, in from_expression
grid=grid, data=expr(**points), label=label, dtype=dtype
File "/home/dperez/.local/lib/python3.8/site-packages/pde/tools/expressions.py", line 487, in call
return self._get_function_cached(single_arg=False)(*args, **kwargs)
File "", line 2, in _lambdifygenerated
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

I use the latest version (pip install git+https://github.com/zwicker-group/py-pde.git

I will appreciate very much your help.

Best,

Daniela

[david-zwicker]

I must admit I have never seen a package from GitHub installed via pip, so I don't know how this works exactly. However, the error message suggests that the fix that we recently implemented is not present in the version of the package that is loaded by default. Could you please double check whether you have the most recent version? in particular, there should be a except ValueError: on line 65 in pde/fields/scalar.py, which fixed the problem.

[DanielaScarpelli]

Dear David,

You were right. I have not installed correctly the new version of the package. Now, the implementation of the Piecewise function is working ok.

I would be very grateful if you could please help me with the following issue: I use the PDE class, as you previously suggested, to define an equation that has coefficients that depends on functions I have previously defined. However, when I run the program I get many errors. I know that the functions I have defined work ok, because when I call them I have no issues. This is part of the code:

from pde import PDE, ScalarField, CartesianGrid
from sympy.abc import x
from sympy import Piecewise
import functions as func

xc = 9770
xd = 10329
xf = 10600

grid = CartesianGrid([[xc, xf]], 50)

print(func.c_tot(xf, 3 * 10 ** (-10)))

f_ini = 1 - ((xf - x) / (xf - xd)) ** 2
p = Piecewise((1 - ((xf - x) / (xf - xd)) ** 2, x >= xd), (0, True))
field = ScalarField.from_expression(grid, f"{p}") # generate initial condition

Define boundary conditions

bc_lower = {"value": 0}
bc_upper = {'type': 'mixed', 'value': 1 / xf, 'const': 0}
#bc = [bc_lower, bc_upper]

Expanded definition of the PDE

yr = 365 * 24 * 60 * 60
c = 3108
#m_accr = 3 * 10(-10)
dif = "c2/(4 * math.pi * func.c_tot(x,3 * 10(-10))) * yr"
mult = " -0.5 * (1/(x**2))"
term_1 = f"({dif}) laplace(s)"
term_2 = f"({dif})*({mult})*s"
eq = PDE({"s": f"{term_1} + {term_2}"}, bc=[bc_lower, bc_upper])
result = eq.solve(field, t_range=10 **3, dt=100)

The messages I get are:
ValueError: Error from parse_expr with transformed code: "(Symbol ('c' )Integer (2 )/(Integer (4 )*Symbol ('math' ).pi *Symbol ('func' ).c_tot (Symbol ('x' ),Integer (3 )*Integer (10 )(-Integer (10 ))))*Symbol ('yr' ))Function ('laplace' )(Symbol ('s' ))+(Symbol ('c' )Integer (2 )/(Integer (4 )*Symbol ('math' ).pi *Symbol ('func' ).c_tot (Symbol ('x' ),Integer (3 )*Integer (10 )(-Integer (10 ))))Symbol ('yr' ))(-Float ('0.5' )(Integer (1 )/(Symbol ('x' )**Integer (2 ))))*Symbol ('s' )"

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
File "/Users/Daniela/PycharmProjects/diffusion_equation/main.py", line 39, in
eq = PDE({"s": f"{term_1} + {term_2}"}, bc=[bc_lower, bc_upper])
File "/Users/Daniela/PycharmProjects/diffusion_equation/venv/lib/python3.7/site-packages/pde/pdes/pde.py", line 135, in init
rhs_expr = ScalarExpression(rhs_str, user_funcs=user_funcs, consts=consts_d)
File "/Users/Daniela/PycharmProjects/diffusion_equation/venv/lib/python3.7/site-packages/pde/tools/expressions.py", line 588, in init
sympy_expr = sympy_parser.parse_expr(expression)
File "/Users/Daniela/PycharmProjects/diffusion_equation/venv/lib/python3.7/site-packages/sympy/parsing/sympy_parser.py", line 1101, in parse_expr
raise e from ValueError(f"Error from parse_expr with transformed code: {code!r}")
File "/Users/Daniela/PycharmProjects/diffusion_equation/venv/lib/python3.7/site-packages/sympy/parsing/sympy_parser.py", line 1092, in parse_expr
rv = eval_expr(code, local_dict, global_dict)
File "/Users/Daniela/PycharmProjects/diffusion_equation/venv/lib/python3.7/site-packages/sympy/parsing/sympy_parser.py", line 908, in eval_expr
code, global_dict, local_dict) # take local objects in preference
File "", line 1, in
AttributeError: 'Symbol' object has no attribute 'pi'

I guess that the syntax to define the equation is wrong, but I do not know how to fix it. I have been dealing for almost a week with this matter.

Many thanks for your time and assistance.

Best,

Daniela

[david-zwicker]

There might be several problems with your code. One certainly is that the mathematical constant pi should be referred by pi instead of math.pi. The string needs to be interpreted by sympy, which has its own syntax. In particular, sympy will not know your function func.c_tot. You need to define it as an actual expression. Alternatively, you could write a custom PDE class, where you can use python functions instead.