Hi, IDK if this is the right place to answer the question, I am solving a 1d wave equation over a bar, if I make the time step size really small, like 10^-5, the max amplitude is 100 times to the time step size of 10^-4 at the same simulation time, both them has the default spatial discretization. Does anyone know what is the issue?
Please post a minimal working example.
Note that depending on the temporal discretization scheme, you Need to fullfill the CFL condition.
For more info see for instance
Hi dokken, since in Fenics, for spatial discretization I only need to write function grad(), if I want to look the code behind the grad(), how to do it? for CFL, I only decrease the time step size, so the CFL will be smaller than 1.
Could you please post a minimal code example illustrating the issue you are having? As you need to discretize both in space and time there are many ways of doing this.
from dolfin import *
import numpy as np
E=3300000000.0
rho=1300.0
c = sqrt(E/rho)
Po=30000000.0
t1=0.00015
t2=0.001
b0=1.0
Length=30.0
# Time variables
dt = 0.00001; t = 0; T = 0.01
mesh = IntervalMesh(1000,0,Length)
V=FunctionSpace(mesh, "Lagrange", 1)
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return near (x[0], 0.0)
class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
return near (x[0], Length)
RB=RightBoundary()
LB=LeftBoundary()
boundaries = MeshFunction("size_t", mesh, 0)
boundaries.set_all(0)
RB.mark(boundaries, 1)
LB.mark(boundaries, 2)
blast=Constant((Po*(1-t/t2)*np.exp(-b0*((t-t1)/t2)))/E)
def blast_func(t):
return (Po*(1-t/t2)*np.exp(-b0*((t-t1)/t2)))/E
ds = Measure("ds", subdomain_data=boundaries)
# 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)
a = u*v*dx + dt*dt*c*c*inner(grad(u), grad(v))*dx
L = 2*u1*v*dx - u0*v*dx + blast*v*ds(2)
bc = DirichletBC(V, 0.0, boundaries, 1)
u=Function(V)
vtkfile1=File('wave/u.pvd')
vtkfile2=File('wave/v.pvd')
while t <= T:
A, b = assemble_system(a, L, bc)
solve(A, u.vector(), b)
velocity = (u - u0) / dt
u0.assign(u1)
u1.assign(u)
t += dt
blast.assign(blast_func(t))
vtkfile1 << (u,t)
vtkfile2 << (project(velocity, V), t)
Both two figures are capture at t=0.01. but with dt =0.00001 and 0.0001
I think I know the reason, because the values of boundary condition at two different timestep are different, but I just didn’t think it will enlarge the result that much.
and another silly question. For total simulation time= 0.01, if I chose the dt=0.001, there will be 10 vtu files, from 000000 to 000009, but if I chose dt=0.0001, there will be 101 vtu files, from 000000 to 000100, why will be like this?
I believe you haven’t scaled your PDE appropriately. Why does your surface integral multiplied by dt**2?
You have a very different boundary condition for each time step if it is note scaled with dt.
Because you are saving a solution per time step, and each solution is saved in a separate file, with a counter corresponding to when it was added. If you only want to have a single file, you need to use XDMFFile
Thank you for your reply dokken, but I don’t quite understand which part you mean for the square of timestep
dt**2
This ds(2) is for my Neumann boundary condition, and this Neumann boundary is a function of time,
When you integrate by parts, you should have a dt**2 times your blast function, as it stems from
oh, not just dt*dt, and c*c
from dolfin import *
import numpy as np
E=3300000000.0
rho=1300.0
c = sqrt(E/rho)
Po=30000000.0
t1=0.00015
t2=0.001
b0=1.0
Length=30.0
# Time variables
dt = 0.0001; t = 0; T = 0.02
mesh = IntervalMesh(1000,0,Length)
V=FunctionSpace(mesh, "Lagrange", 1)
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return near (x[0], 0.0)
class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
return near (x[0], Length)
RB=RightBoundary()
LB=LeftBoundary()
boundaries = MeshFunction("size_t", mesh, 0)
boundaries.set_all(0)
RB.mark(boundaries, 1)
LB.mark(boundaries, 2)
blast=Constant((Po*(1-t/t2)*np.exp(-b0*((t-t1)/t2)))/E)
def blast_func(t):
return (Po*(1-t/t2)*np.exp(-b0*((t-t1)/t2)))/E
ds = Measure("ds", subdomain_data=boundaries)
# 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)
a = u*v*dx + dt*dt*c*c*inner(grad(u), grad(v))*dx
L = 2*u1*v*dx - u0*v*dx + dt*dt*c*c*blast*v*ds(2)
bc = DirichletBC(V, 0.0, boundaries, 1)
u=Function(V)
vtkfile1=File('wave/u.pvd')
vtkfile2=File('wave/v.pvd')
while t <= T:
A, b = assemble_system(a, L, bc)
solve(A, u.vector(), b)
velocity = (u - u0) / dt
u0.assign(u1)
u1.assign(u)
t += dt
blast.assign(blast_func(t))
vtkfile1 << (u,t)
vtkfile2 << (project(velocity, V), t)
Both two figures are capture at t=0.02, but with dt =0.00001 and 0.0001
Is this an acceptable error?
Sorry, I didn’t make myself clear, I mean if dt=0.0001, it should be 100 files which are from 000000 to 000099, just like when dt=0.001. But now it is from 00000 to 000100, why will be like this?
Please format your code appropriately.
I believe you should check the solution for a range of time steps, to observe how you solution is dependening on the temporal discretization.
Depending on what your boundary condition is (you have not stated your strong form), the boundary term should not be scaled by c*c.
Please carefully formulate your strong problem and derive your weak form.
The PDE is
\frac{\partial^2u}{\partial t^2}=c^2\frac{\partial^2u}{\partial x^2}
ICs are
u(x,0)=0
u_t(x,0)=0
BCs are
u_x(0,t)=f(t)
u(L,t)=0
where x\in[0, L]
Second order temporal dicretization is
(\frac{\partial^2u}{\partial t^2})^{n+1}=\frac{u^{n+1}-{2u^n}+u^{n-1}}{\Delta t^2}
so we have
u^{n+1}-{2u^n}+u^{n-1}=\Delta t^2 c^2(\frac{\partial^2u}{\partial x^2})^{n+1}
then the variational form is
u^{n+1}v-\Delta t^2 c^2(\frac{\partial^2u}{\partial x^2})^{n+1}v=2u^n v-u^{n-1}v
then apply the boundary condition by integration by parts to the (\frac{\partial^2u}{\partial x^2})^{n+1}v, so there will be a scale factor \Delta t^2 c^2 to the Neuman boundary
Continue refining you temporal discretization, and you will observe convergence of your solution.
Any suggestion for improving the eps? this is the print for time,
