Enforce Dirichlet Boundary Condition manually on solution function

1

Hey,

this is probably super easy, but I can’t seem to find a suitable solution. I implemented my own Newton sovler. My boundary conditions of my system change with time. In every iteration in my newton solver I want to reuse the solution from my previous iteration, however since in my Newton solver the boundary conditions are only applied as homogenous ones, their values are not updated with time. Therefore I want to set the values at the boundaries manually (update u_k). How is this easiest done and still applicable for generic meshes?
This is a MWE hopefully showing my problem:

from fenics import *
def CustomNewton(F, u,bc, V):
   
    J = derivative(F,u)
    U_inc  = Function(V)  

    #homogenize boundary conditions 
    bc.homogenize()

    # initialize the variables for the loop
    nIter     = 0
    eps       = 1e14

    # Initialize the solver for the linear system
    solver = KrylovSolver('gmres', 'ilu')
    while eps > 1e-15 and nIter < 20:        
        # Newton iterations
        nIter += 1

        A = assemble(J, keep_diagonal=True)
        b = assemble(-F)

        bc.apply(A,b)

        A.ident_zeros() # since we use mixedElements which are not all defined on the whole domain

        eps = norm(U_inc)
        solver.solve(A, U_inc.vector(), b)     # Determine step direction and size
        
        u.assign(u+0.5*U_inc)
        print ('numIter = ',nIter,', eps = ' ,eps)

mesh = UnitSquareMesh(64,64)
V    = FunctionSpace(mesh, 'Lagrange', 1)

u_k   = Function(V)
u_k_1 = Function(V)
v = TestFunction(V)
dt = 0.1
        
F = (u_k-u_k_1)/dt*v*dx - dot(u_k*grad(u_k), grad(v))*dx - 4*u_k*v*dx

bc_val = Expression('t+pow(x[0],2)+pow(x[1],2)', t=0, degree=2)
bc = DirichletBC(V,bc_val , 'on_boundary')

u_k.interpolate(Expression('pow(x[0],2)+pow(x[1],2)', degree=0))   # load initial solution
u_k_1.interpolate(Expression('pow(x[0],2)+pow(x[1],2)', degree=0)) 

for k in range(0,10):
    t = k*dt
    CustomNewton(F, u_k,bc, V)
    # test value at (1,1)
    print('Value at (1,1) ' + str(u_k(Point(1.,1.))) + ' (it should be ' + str(2+t) + ')')
    #refresh Boundary
    bc_val.t = t
    u_k_1.assign(u_k)
    # Here: update the boundary values of u_k in order to simulate the dynamic system

Thank you for your time :slight_smile:
Best Regards

Emanuel

2

Ok, I found it: It’s as easy as I thought it was. Using

bc.apply(u_k.vector())

one can enforce the boundaries onto u_k. You then just have to create a homogenized version of the bc and a non-homogenized version.
For reference the working implementation

from fenics import *

def CustomNewton(F, u,bc, V):
   
    J = derivative(F,u)
    U_inc  = Function(V)  

    #homogenize boundary conditions 
    bc.homogenize()

    # initialize the variables for the loop
    nIter     = 0
    eps       = 1e14

    # Initialize the solver for the linear system
    solver = KrylovSolver('gmres', 'ilu')
    while eps > 1e-15 and nIter < 20:        
        # Newton iterations
        nIter += 1

        A = assemble(J, keep_diagonal=True)
        b = assemble(-F)

        bc.apply(A,b)

        A.ident_zeros() # since we use mixedElements which are not all defined on the whole domain

        eps = norm(U_inc)
        solver.solve(A, U_inc.vector(), b)     # Determine step direction and size
        
        u.assign(u+0.5*U_inc)
        print ('numIter = ',nIter,', eps = ' ,eps)

mesh = UnitSquareMesh(64,64)
V    = FunctionSpace(mesh, 'Lagrange', 1)

u_k   = Function(V)
u_k_1 = Function(V)
v = TestFunction(V)
dt = 0.1
        
F = (u_k-u_k_1)/dt*v*dx - dot(u_k*grad(u_k), grad(v))*dx - 4*u_k*v*dx

bc_val = Expression('t+pow(x[0],2)+pow(x[1],2)', t=0, degree=2)
bc = DirichletBC(V,bc_val , 'on_boundary')
bc_Hom = DirichletBC(V,bc_val , 'on_boundary')

u_k.interpolate(Expression('pow(x[0],2)+pow(x[1],2)', degree=0))   # load initial solution
u_k_1.interpolate(Expression('pow(x[0],2)+pow(x[1],2)', degree=0)) 

for k in range(0,10):
    t = k*dt
    #refresh Boundary
    bc_val.t = t
    bc.apply(u_k.vector())
    CustomNewton(F, u_k,bc_Hom, V)
    # test value at (1,1)
    print('Value at (1,1) ' + str(u_k(Point(1.,1.))) + ' (it should be ' + str(2+t) + ')')
    
    
    u_k_1.assign(u_k)