Consider a very simplified 1d wave equation, and with two easy initial conditions
u(x,0)=0,
u_t(x,0)=0,
Well assuming you have a time discretization then you need to define three functions u and u_n, u_nn.
u will be used for time n+2 and u_n for time n and u_nn for n+1. These are to be able to discretize the second order derivative you have in the wave equation.
In Fenics given a proper vector space V the terms
u = Function(V)
u_n = Function(V)
u_n n= Function(V)
initialize them with zero values. As for the second question, you can for example decrease your time steps. Or take a look at this.
Thank you nami, I have another question, for this 1d wave equation, I also have two mixed boundary conditions, namely
u(L,t)=0, where x belongs to 0 to L,
u_x(0,t)=f(t), f(t) is a function of t, like exp(t)
How to handle these boundary conditions?
BCs like these have to be incorporated in the weak formulation. Namely
\int_{\Omega} u_{xx} v dX = \int_{\partial \Omega} \frac{\partial u}{\partial \hat{n}} v dS - \int_{\Omega} u_{x} v_x dX
where \hat{n} is the unit outward normal vector at the boundaries (in your 1D case it is 1 or -1 depending on whether you are at L or 0). The function v is a test function and \Omega is your domain. The right hand side comes form a general integration by parts and the boundary integral is where you should apply the BCs.
Yes, I got this.
I mean is there an existed 1d mixed boundary conditions code (similar), so I can easily modify it? I really don’t get the syntax of Fenics
Oh sorry. So you are looking for a source term on one side. Then something like this should work:
from dolfin import *
import numpy as np
c =1000
mesh = IntervalMesh(100,0,1.0)
V=FunctionSpace(mesh, "Lagrange", 1)
class Right(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[0],1)
boundaries = MeshFunction("size_t", mesh,0)
boundaries.set_all(0)
# Now mark your Neumann boundary
rightbound = Right()
rightbound.mark(boundaries, 1)
ds=Measure('ds',subdomain_data=boundaries)
# Time variables
dt = 0.001; t = 0; T = 1
# Previous and current solution
u1= interpolate(Constant(0.0), V)
u0= interpolate(Constant(0.0), V)
# Variational problem at each time
u = TrialFunction(V)
v = TestFunction(V)
left = CompiledSubDomain("(std::abs(x[0] - 0.0) < DOLFIN_EPS) && on_boundary")
a = u*v*dx + dt*dt*c*c*inner(grad(u), grad(v))*dx
L = 2*u1*v*dx-u0*v*dx+ exp(t)*v*ds
bc = DirichletBC(V, 0, left)
A, b = assemble_system(a, L, bc)
u=Function(V)
vtkfile=File('wave/u.pvd')
while t <= T:
A, b = assemble_system(a, L, bc)
solve(A, u.vector(), b)
u0.assign(u1)
u1.assign(u)
t += dt
vtkfile<<(u,t)
I got it from here and changed it a bit to match your question.
Thank you so much nami, and I have a stupid question when I use paraview to view the results of this code, it has this
That is strange. I don’t have that problem. Also I edited the code a second time. Make sure you are using my last edit. There were a couple of issues.
What version of ParaView you are using?
I am using 5.4.1 64-bit for windows. But I have an older version on my Linux computer and it works there as well.
Oh I see,
Hi, if I also want to solve for velocity, which is u_t(x,t), what is the simplest way to implement it based on this code?
Do you mean you want to know what your velocity is based on the displacements you get our you want to break the problem into a system that involves velocity? I feel like you mean the former so you just use a discrete formula right before u0.assign(u1), like:
velocity = (u-u0)/dt
First, you need to create a second file, if you do not want to overwrite your velocity on your already saved displacements. However, this is not why you are getting the error and the error is coming from the division. As you can see in the error. To fix that, first project the velocity and then save it:
vtkfile << (project(velocity,V),t)Thank you nami, if I may ask, why I need to project the velocity, I think it should be exactly as same as displacement.
u and u0 are functions so Fenics already knows the values of them at each node. However, once you mix it with a non-function object such as dt it isn o longer the case. So the projection defines the values of the new quantity at each node.
Thank you, and for the previous code, x in [0, L]
the D BC is u(L,t)=0, so it should be the right side, I try to modify the code but get an error,
bc = DirichletBC(V, 0, Right)
A, b = assemble_system(a, L, bc)
Why here cannot be
(u-u1)/dt
u1 is closer than u0 to u, but I test, this one, it is all zero