Time depended Neumann boundary

Miraboreasu · February 17, 2021, 4:00pm

I have a time depended Neumann boundary, exp(t) for 1d wave equation, I code it into the weak form, but when I put it into the while loop, it got an error

from dolfin import *
import numpy as np

c=1000
Length=100.0

mesh = IntervalMesh(1000,0,Length)
V=FunctionSpace(mesh, "Lagrange", 1)

class Right(SubDomain):
    def inside(self,x,on_boundary):
        return on_boundary and near(x[0],Length)

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.0001; t = 0; T = 0.02

# 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)

bc = DirichletBC(V, 0, rightbound)


u=Function(V)
vtkfile1=File('wave/u.pvd')
vtkfile2=File('wave/v.pvd')
while t <= T:

    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
  

    A, b = assemble_system(a, L, bc)
    solve(A, u.vector(), b)

dokken · February 17, 2021, 4:03pm

Please follow the guidelines (Read before posting: How do I get my question answered?)
and post a complete minimal code example and include the whole error message.

Miraboreasu · February 17, 2021, 4:10pm

Sorry, I didn’t see the rest part of the error. The post has been modified.

dokken · February 17, 2021, 4:16pm

You are redefining the trial function u before you are defining your variational form, and Therefore you obtain an error message. As a general rule, you should try to define the variational formulation outside for loops.

Miraboreasu · February 17, 2021, 4:17pm

If I put variational formulation out the while loop, how to update the time for the boundary?

dokken · February 17, 2021, 4:23pm

Define either an Expression or a Constant which can you can update the value of inside the for-loop. See for instance: Time-integration of elastodynamics equation — DOLFIN documentation
or Bitbucket

Miraboreasu · February 17, 2021, 7:03pm

I am follow the How to apply time varying boundary conditions coming from a function?

I test this with two def, one is exp(t), another one is a long equation, but except t, all other letters are defined, which are numbers.

blast=Constant((Po*(1-t/t1)*exp(-b*((t-t2)/t1)))/E)
def blast_func(t):
    return (Po*(1-t/t1)*exp(-b*((t-t2)/t1)))/E

pexp=Constant(exp(t))
def complicate_func(t):
    return exp(t)

a = u*v*dx + dt*dt*c*c*inner(grad(u), grad(v))*dx
L = 2*u1*v*dx - u0*v*dx + pexp*v*ds

bc = DirichletBC(V, 0, rightbound)
#A, b = assemble_system(a, L, bc)

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
    pexp.assign(complicate_func(t))
    blast.assign(blast_func(t))
    vtkfile1<<(u,t)
    vtkfile2 << (project(velocity, V), t)

dokken · February 17, 2021, 7:13pm

You need to make sure that the return type of blast_func is a float. Currently you are using ufl.exp (as it is import with dolfin). You need to use numpy.exp to make the return type a float.

For further questions, please create a minimal code that reproduces the error (and include no lines that are not required to reproduce the error). The code should be possible to copy paste into an empty file and run directly.

Miraboreasu · February 17, 2021, 7:27pm

I replace the exp with np.exp

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
b=1.0
Length=100.0

mesh = IntervalMesh(1000,0,Length)
V=FunctionSpace(mesh, "Lagrange", 1)

class Right(SubDomain):
    def inside(self,x,on_boundary):
        return on_boundary and near(x[0],Length)

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.0001; t = 0; T = 0.02

# 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)

blast=Constant((Po*(1-t/t1)*np.exp(-b*((t-t2)/t1)))/E)
def blast_func(t):
    return (Po*(1-t/t1)*np.exp(-b*((t-t2)/t1)))/E

pexp=Constant(exp(t))
def complicate_func(t):
    return exp(t)

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

bc = DirichletBC(V, 0, rightbound)
#A, b = assemble_system(a, L, bc)

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
    #pexp.assign(complicate_func(t))
    blast.assign(blast_func(t))
    vtkfile1 << (u,t)
    vtkfile2 << (project(velocity, V), t)

dokken · February 17, 2021, 7:32pm

Again, you keep on redeclaring variables. In this case it is b

that is redeclared as

Call the first b something different, for instance b_.

Miraboreasu · February 18, 2021, 6:39am

Thank you dokken, another boundaries question, I learn another way to set the mixed boundary conditions(22. Poisson equation with multiple subdomains — FEniCS Project), namely

for x in [0,Length]

u=(Length, t)=0
u_x(0, t)=func(t)

The first (previous) version is

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.1

mesh = IntervalMesh(1000,0,Length)
V=FunctionSpace(mesh, "Lagrange", 1)

class Right(SubDomain):
    def inside(self,x,on_boundary):
        return on_boundary and near(x[0], Length)

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)


# 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)


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

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

bc = DirichletBC(V, 0, rightbound)


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)

The second version is

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.1

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

# 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)

In the first one, the wave will rebound from the right side, the second is all zero, I don’t know where the problem is in my second version code.

dokken · February 18, 2021, 7:09pm

You seem to be not reading through the tutorial carefully (and note that you are looking at a very old tutorial, the latest version can be found here:Poisson equation with multiple subdomains — DOLFIN documentation).
You need to define the integration measure that is related to your boundary marker:

ds = Measure("ds", subdomain_data=boundaries)

Please read carefully through the tutorial before asking questions.

Topic		Replies	Views
Setting neumann boundary condition for mixed function spaces Errors	4	588	June 1, 2022
Solving 1 space 1 time equation with Neumann bc vs Dirichlet bc	1	563	May 2, 2020
Cannot get 1D time dependent solution General	1	363	November 9, 2022
Time dependent problem with weak form depending on the solution of a secondary problem variational formulation	0	280	October 12, 2022
Evaluating time dependent Neumann boundary condition using tensors I/O	5	1093	January 16, 2022

Time depended Neumann boundary

Related topics