Updating temporal boundary conditions

variational formulation
1

Hello, i’m solving this Navier-Stokes-like equations but when i plot the solution (with Visit) i see that the boundary conditions don’t update like they should (the arrows may turn because it’s a cosine, sin function and they don’t do that), can you find the problem? Here’s the code:

%matplotlib inline
import numpy as np
import math
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
f=0;
g=9.81;

c_d=Constant(0.0);
α_h=Constant(0.0);

z0=1.0;
ze=1.0;
L=1.0;
om1=f/2+sqrt(f**2/4+2*g*z0/L**2);

num_steps=1000;#pasos de tiempo
T=2*pi/om1;
dt=T/num_steps;
k  = Constant(dt)
mod=10 

center = Point(0.0,0.0)
horizontal_semi_axis = L
vertical_semi_axis = L
domain = Ellipse(center,
        horizontal_semi_axis,
        vertical_semi_axis)
resolution = 5
mesh = generate_mesh(domain,resolution)

ufile = File("SWE_results/velocity.pvd")
ηfile = File("SWE_results/elevation.pvd")

V = VectorFunctionSpace(mesh, "Lagrange", 1)
Q = FunctionSpace(mesh, "Lagrange", 1)

u = TrialFunction(V)
η = TrialFunction(Q)
v = TestFunction(V)
q = TestFunction(Q)

u_n = Function(V)
u1  = Function(V)
η_n = Function(Q)
η1  = Function(Q)

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
class boundary(SubDomain):
def inside(self, x, on_boundary):
    return on_boundary
boundary = boundary()
boundary.mark(boundary_markers,0)
t=0.0
η_D=Expression('2*ze*z0/L*(x[0]/L*cos(om1*t)-x[1]/L*sin(om1*t)-ze/2/L)', degree=2, ze=ze, z0=z0, L=L, om1=om1, t=t)
u_D=Expression(('-ze*om1*sin(om1*t)','-ze*om1*cos(om1*t)'), degree=2, ze=ze, om1=om1, t=t)

bcη = DirichletBC(Q, η_D, boundary)
bcu = DirichletBC(V, u_D, boundary)

bcη=[bcη]
bcu=[bcu]

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
n = FacetNormal(mesh)

h=project(Expression('z0*(1.-pow(x[0],2)/pow(L,2)-pow(x[1],2)/pow(L,2))', degree=2, z0=z0, L=L),Q)

F1 = inner((u - u_n), v)*dx + k*inner(α_h*grad(u_n),grad(v))*dx + k*inner(grad(u_n)*u_n, v)*dx + k*k/4*inner(grad(u_n)*u_n,grad(v)*u_n)*dx
a1 = lhs(F1)
L1 = rhs(F1)
a2 = inner(η,q)*dx + k*k/4*g*inner((η_n+h)*grad(η),grad(q))*dx
L2 = inner(η_n,q)*dx - k*k/4*g*inner((η_n+h)*grad(η_n),grad(q))*dx - k/2*inner(div((η_n+h)*u1),q)*dx - k/2*inner(div((η_n+h)*u_n),q)*dx
a3 = inner(u,v)*dx
L3 = inner(u1,v)*dx - k/2*g*inner(grad(η1) - grad(η_n),v)*dx

prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"
parameters['krylov_solver']['nonzero_initial_guess'] = True

ufile = File("results/velocity.pvd")
ηfile = File("results/elevation.pvd")

η_i=project(Expression('2*ze*z0/L*(x[0]/L*cos(om1*t)-x[1]/L*sin(om1*t)-ze/2/L)', degree=2, ze=ze, z0=z0, L=L, om1=om1, t=0.0),Q)
u_i=project(Expression(('-ze*om1*sin(om1*t)','-ze*om1*cos(om1*t)'), degree=2, ze=ze, om1=om1, t=0.0),V)

u_n.assign(u_i)
η_n.assign(η_i)
t = 0.0
for m in range(num_steps):
    η_i.t = t    
    A1 = assemble(a1)
    [bc.apply(A1) for bc in bcu]    
    begin("Computing tentative velocity")
    b1 = assemble(L1)
    [bc.apply(A1, b1) for bc in bcu]
    solve(A1, u1.vector(), b1, "bicgstab", "default")
    end()

    A2 = assemble(a2)    
    [bc.apply(A2) for bc in bcη]
    begin("Computing elevation correction")
    b2 = assemble(L2)
    [bc.apply(A2, b2) for bc in bcη]
    [bc.apply(η1.vector()) for bc in bcη]
    solve(A2, η1.vector(), b2, "bicgstab", prec)
    end()

    A3 = assemble(a3)    
    begin("Computing velocity correction")
    b3 = assemble(L3)
    [bc.apply(A3, b3) for bc in bcu]
    solve(A3, u1.vector(), b3, "bicgstab", "default")
    end()

    if (m%mod==0): 
          ufile << u1
          ηfile << η1
    u_n.assign(u1)
    η_n.assign(η1)
    t += dt
2

I don’t see anywhere in your time loop u_D.t=t. That’s probably the issue.

3

Well that’s right, thank you!