[solve command syntax error] 2D Wave equation

variational formulation
1

Hello, I am trying to simulate a 2D wave equation on a sample square mesh but when I assemble the lhs, rhs or use the solve command I get syntax error:

Blockquote
File “”, line 11
solve(A==K, u_s, solver_parameters=dict(linear_solver=‘gmres’, preconditioner=‘ilu’))
^
SyntaxError: invalid syntax

I have 2 trial function (u_1 is the pressure and u_2 is the velocity vector field). I have 2 equations (first discretized with implicit euler and the second with explicit euler), which I summed them up. This is what I am trying to implement in FEniCS:

\text{(Continuous problem): } \begin{array}{ll} \chi_s\frac{\partial}{\partial t}\int_\Omega \bar v_1 \cdot \bar u_1\,d\Omega + \mu_0 \frac{\partial}{\partial t}\int_\Omega \bar v_2 \cdot \bar u_2\,d\Omega = \int_\Omega \text{grad}\bar v_1 \cdot \bar u_2\,d\Omega + \int_{\partial \Omega}\bar v_1\cdot \bar u_2\cdot \boldsymbol n\, dS -\int_\Omega \bar v_2\cdot \text{grad}\bar u_1\,d\Omega \end{array}
\text{(Time-discretized): } \begin{array}{ll} \chi_s\int_\Omega \bar v_1 \cdot \frac{u_1^{n+1}-u_1^n}{\Delta t}\,d\Omega + \mu_0\int_\Omega \bar v_2 \cdot \frac{u_2^{n+1}-u_2^n}{\Delta t}\,d\Omega = \int_\Omega \text{grad}\bar v_1 \cdot \bar u_2^{n+1}\,d\Omega + \int_{\partial \Omega}\bar v_1\cdot \bar u_2^{n+1}\cdot \boldsymbol n\, dS -\int_\Omega \bar v_2\cdot \text{grad}\bar u_1^n\,d\Omega \end{array}

I tried different implementations but also assigning lefthand side and righthand side return me the same error of syntax.

F = v1*u1*dx - v1*u_n.split()[0]*dx + dot(v2, u2)*dx - dot(v2, u_n.split()[1]*dx -dt*dot(grad(v1), u2)*dx +dt*dot(v2, grad(u_n.split()[0]))*dx + dt*v1*dot(u_b.split()[1], n)*ds

a, L = lhs(F), rhs(F)
solve(a == L, u_s, bc)

File “”, line 15
a, L = lhs(F), rhs(F)
^
SyntaxError: invalid syntax

Here is my code:

from fenics import *
import numpy as np

# Constants 
dt = 0.0005
FT = 0.39 # final time
t = 0

# Mesh
N = 20
M = 20
mesh = RectangleMesh(Point(0,0), Point(1,1), N, M, "left")
n = FacetNormal(mesh)

# Finite element spaces
P1 = FiniteElement("P-", triangle, 1, 0)
RT1 = FiniteElement("RTE", triangle, 1)
Lmix = FunctionSpace(mesh, P1*RT1)

# Trial functions
u_1, u_2 = TrialFunctions(Lmix)

# Test functions
v_1, v_2 = TestFunctions(Lmix)

# Functions - previous step - boundary inputs
u_n = Function(Lmix) # previous step (respect to n+1)
u_b = Function(Lmix) # boundary inputs
u_s = Function(Lmix) # solution (to be updated at each time step in the iteration loop)

# Forms
A = v_1*u_1*dx + dot(v_2, u_2)*dx - dt*dot(grad(v_1), u_2)*dx

# Compute the solution at each time step

while(t < 0.05):
    # Update boundary conditions
    g = Expression(("0", "sin(pi*t)*near(x[0],0)*(t<1)", "sin(pi*t)*near(x[1],1)*(t<1)"), degree = 2,t = t)
    u_b = project(g, Lmix)
    
    K = v_1*u_n.split()[0]*dx + dot(v_2, u_n.split()[1]*dx - dt*dot(v_2, grad(u_n.split()[0]))*dx - dt*v_1*dot(u_b.split()[1], n)*ds

    #Solve linear system
    solve(A==K, u_s, solver_parameters=dict(linear_solver='gmres', preconditioner='ilu'))
    
    
    # Assign solution to previous step vector
    u_n.assign(u_s)
    print(t)
    
    # Plot u_1
    c = plot(u_s.split()[0])
    plt.colorbar(c)
    plt.show()
    
    # Update time
    t += dt

Sorry for my bad explanation, it ismy first post and I am beginning to learn how to use the forum. If you need more information, please ask and I will give more details.
Any help or suggestion will be very appreciated.
Thank you very much!

Niccolò

2

You forgot to close the parenthesis of the first dot() in the definition of K.
Also it is not necessary to redefine the expression g and K on every time iteration.

1 Like
3

To complete what dajuno is saying, to update your BC at each timestep just do g.t = t

1 Like
4

@dajuno @Thibaut
Thank you very much for your observations and support!
The syntax mistake of the brackets was really the point that broke the code.

Concerning your idea to define boundary and variational forms outside the loop (did I understand correctly?), I am trying to implement it and, also, I changed a bit the notation for the boundary conditions definition.
I did not get any error but in the simulation the system does not evolve at all, everything is kept at zero, so I think I did not understand completely your point.

Maybe, should I define the initial value for u_n (solution and time n) and u_b (boundary input) also in the loop, or even bc?

Here is the “new” code:

# Trial functions
u1, u2 = TrialFunctions(Lmix)

# Test functions
v1, v2 = TestFunctions(Lmix)

# Functions - previous step - boundary inputs
u_n = Function(Lmix) # previous step (respect to n+1)
u_b = Function(Lmix) # boundary inputs
u_s = Function(Lmix) # solution (to be updated at each time step in the iteration loop)

# Define boundaries
def boundary(x):
        return near(x[0],0) or near(x[0],1) or near(x[1],0) or near(x[1],1)
    
#Boundary conditions
u_D = Expression(("0", "sin(pi*t)*near(x[0],0)*(t<1)", "sin(pi*t)*near(x[1],1)*(t<1)"), degree = 2,t = t)
bc = DirichletBC(Lmix, u_D, boundary)

# Define initial value
u_n = interpolate(u_D, Lmix)
u_b = project(u_D, Lmix)

# Variational problem
F = v1*u1*dx - v1*u_n.split()[0]*dx + dot(v2, u2)*dx - dot(v2, u_n.split()[1])*dx - dt*dot(grad(v1), u2)*dx + \
	dt*dot(v2, grad(u_n.split()[0]))*dx + dt*v1*dot(u_b.split()[1], n)*ds

# Compute the solution at each time step

while(t < FT):
  
    # Update Boundary conditions
    u_D.t = t
      
    a, L = lhs(F), rhs(F)
    solve(a == L, u_s, bc)
                                                              
    # Assign solution to previous step vector
    u_n.assign(u_s)
    print(t)
    
    # Plot u_1 (pressure)
    c = plot(u_s.split()[0])
    plt.colorbar(c)
    plt.show()
    
    # Update current time
    t += dt
5

Ok I quickly checked your formulation. I don’t know the problem and cannot really judge the setting. However, a couple of things. The sign of dt*v1*dot(u_b.split()[1], n)*ds is wrong, right? Why the u_b term anyways? Doesn’t the DirichletBC handle this? Judging from your formulation, this should rather be u2?
If you need to work with this function u_b, and it’s supposed to equal the boundary Expression u_D, you have to project it everytime you update u_D.t. This is missing in your code.
Hope this helps

6

Hello, sorry for the delay!

  1. You’re right! The sign was wrong.

  2. Yes, from the formulation it should be u2 rather than ub, but since that term in evaluated on the boundary domain, that will be our “control input”. To be honest I am working on a project that was begun from another guy and I kept his notation (he introduced u_b and I am not sure if it is correct anyway). I am not sure also about the mechanism of DirichletBC, because if (as you suggest) they already take it into account, then maybe this is not a good solution; but the code seemed to work so I tried to change as little as possible.

  3. I tried your idea but it did not update anything (really I don’t know why), so for the moment I probably will keep the code as i was at the beginning.

Now I am trying to transpose the problem into 3D, and I am having a hard time implementing and updating these BCs. So I will open another topic :sweat_smile:

For now, thank you again for your help!

7

Just a personal comment, I guess without really understanding what’s going on in that code and with the boundary conditions, it doesn’t make much sense turning to 3D :wink:
On the contrary, I think it’d be a good idea to simplify the problem as much as possible to eliminate possible causes of the errors and isolate the issue.