How to apply Dirichlet BC to some given degree of freedom


112
views
0
9 weeks ago by
Hi everyone,
I am solving a PDE problem on the domain Omega = (0,Lx) x (0, Ly).
-div (h^2 * grad(u)) = f in Omega.
                            u = 0 at x = 0 or x = 1.
suppose that h is given.
We consider Dirichlet boundary condition on the left, x=0, and the right, x = 1. We set the upper bound for this solution, say u <= ub = 0.2 where h = 1 and f = 2. After solving this PDE problem, there are some nodes of the mesh reaching the upper bound. Let us call the set UAS, which contains some degree of freedom of the problem.
Now I would like to solve the self-adjoint problem PDE problem i.e. consider the PDE with the same settings and  an additional Dirichlet boundary condition on the set UAS, u = 0 on UAS.
My question is that how can we apply the Dirichlet condition on some given  degree of freedom (elements of UAS)? Note that I use IndexSet to mark the UAS set.

Some lines of code that I am doing are attached. Thank you so much.

Best regards,
Nha
"""
    FEniCS program:
   
    -div(q(h,p)*grad(u)) = f   in (0, Lx) x (0, Ly)
                       u = uD on the boundary.
    
    q(h,p) = h^p
    i.e. p = 2 and mu =2
    """
# This programe is a modification from the code of Dr. Maurini to study a PDE problem.

# Follow the links for more infomation
# http://www.lmm.jussieu.fr/~corrado/_site/codes/
# https://bitbucket.org/cmaurini/varfrac_for_cism
#==========================================================
#==========================================================
# Copyright (C) 2012-2013 Corrado Maurini.
# Licensed under the GNU LGPL Version 3.
#
#==========================================================
#                      Import libraries
#==========================================================
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

from dolfin import *
import numpy as np
from numpy import linalg as LA
from petsc4py import PETSc
try:
    from petsc4py import PETSc
except:
    print ("\n*** Warning: you need to have petsc4py installed for this module to run\n")
    print ("\nExiting.\n")
    exit()

if not has_petsc4py():
    print ("\n*** Warning: Dolfin is not compiled with petsc4py support\n")
    print ("\nExiting.\n")
    exit()

#==========================================================
#==========================================================
#                   The main function
#==========================================================
#==========================================================
def main():
    parameters["linear_algebra_backend"] = "PETSc"
    #==========================================================
    #                        Parameters
    #==========================================================

    # Some constants
    mu = 2.; p = 2.; tol = 1e-14; zero = Constant(0.0);

    ubVal = 0.2; lbVal = -100;

    # Create mesh and define function space
    Lx = 1.; Ly = 1.; Nx = 4; Ny = 4;
    # The mesh
    mesh = RectangleMesh(Point(0,0), Point(Lx, Ly), Nx, Ny, "right/left")

    # Degree of the varialbes
    u_degree = 2; grad_u_degree = 1; h_degree = 1;
    # Define functional spaces
    Vh = FunctionSpace(mesh, 'P', h_degree)
    Vu = FunctionSpace(mesh, 'P', u_degree)
    Vu_grad = FunctionSpace(mesh, 'P', grad_u_degree)


    # Define Dirichlet boundary conditions Rec
    def BoundaryX0(x, on_boundary):
        return on_boundary and near(x[0], 0., tol)
    
    def BoundaryX1(x, on_boundary):
        return on_boundary and near(x[0], 1., tol)
    # Apply B.C to Vu
    bcX0 = DirichletBC(Vu, zero, BoundaryX0)
    bcX1 = DirichletBC(Vu, zero, BoundaryX1)
    bcs = [bcX0, bcX1]

    # Define variational problem
    u = TrialFunction(Vu)
    usol = Function(Vu)
    v = TestFunction(Vu)
    h = Function(Vh)

    # The RHS of the PDE
    fVal = 2.
    f = Expression('fVal', degree=1, fVal = fVal)

    # Initial guess for h
    h0Val = 1.
    h0 = Expression('h0Val', degree=1, h0Val = h0Val)
    h.assign(interpolate(h0, Vh))

    a = (h**p)*inner(grad(u), grad(v))*dx
    L = f*v*dx

    # Define the upper and lower bounds
    ub = interpolate(Constant(ubVal), Vu)
    lb = interpolate(Constant(lbVal), Vu)
    #-----------------------------
    # Solution with TAO (linear system formulation, working only for quadratic energies)
    #-----------------------------
    # Assemble the linear system
    A, b = assemble_system(a, L, bcs)
    # redefine the functional space for the displacement field, u
    u = Function(Vu)

    solver = TAOLinearBoundSolver("gpcg","gmres")
    # Set some parameters
    solver.parameters["monitor_convergence"]=True
    solver.parameters["report"]=True


    # Solve the problem
    solver.solve(A, u.vector(), b , lb.vector(), ub.vector())
    info(solver.parameters,True)

    ##==========================================================
    ##                  Find the Active set: Begin
    ##==========================================================
    u_array = u.vector().get_local()
    ub_array = ub.vector().get_local()

    UAS_array = np.where(u_array[:] >(ub_array[:] - 1e-5))
    UAS = PETSc.IS().createGeneral(UAS_array)
    UAS.view()
    
#    ##==========================================================
#    ##                  Find the Active set: End
#    ##==========================================================

    # The maximum of u
    u_max = np.abs(u.vector().get_local()).max()
    print('\n u_max = \n', u_max)

    h_max = np.abs(h.vector().get_local()).max()
    print('\n h_max = \n', h_max)

    # plot the solution
    plot(u, title= "Displacement")
    plot(mesh)
    vtkfile_u = File('Output/displacement.pvd')
    vtkfile_u << u

if __name__ == '__main__':
    main()
​
Community: FEniCS Project

1 Answer


3
9 weeks ago by
There may be a smarter way to do this through DOLFIN, but one solution is to just apply the BCs yourself with petsc4py.  Given the desired homogeneous Dirichlet BC DoFs as an index set, you would modify the LHS like
as_backend_type(A).mat().zeroRowsColumns(bc_dofs)
as_backend_type(A).mat().assemblyBegin()
as_backend_type(A).mat().assemblyEnd()​
where A is the LHS matrix and bc_dofs is an IS.  (Note that, despite the name of zeroRowsColumns(), its default behavior is to set the diagonal entry to 1, as desired.)  Then you can modify the RHS like
as_backend_type(b).vec().setValues(bc_dofs,zeros(bc_dofs.getLocalSize()))
as_backend_type(b).vec().assemblyBegin()
as_backend_type(b).vec().assemblyEnd()
where b is the RHS vector, and the zeros() function is imported from numpy.  
Thank you so much for your help, Kamensky. I appreciate it.
I thought about that, however, I could not able to express my thoughts in terms of Fenics's language since Fenics is still new to me.
However, I did try another approach and it does work for now. What i did is to set upper bound and lower bound at these degrees of freedom to be zero. The results turn out to be good.
Thank you so much for your help.
Best wishes,
Nha
"""
    FEniCS program: Optimization with PDE constraint
    
    Min_h {J(h) = \int_\Omega f u dx + \mu * \int_\Omega h^p dx}
    s.t
    -div(q(h,p)*grad(u)) = f   in the unit square.
    u = uD on the boundary.
    
    q(h,p) = h^p
    i.e. p = 2 and mu =2
    """
# This programe is a modification from the code of Dr. Maurini to study a PDE constrained optimization problem.

# Follow the links for more infomation
# http://www.lmm.jussieu.fr/~corrado/_site/codes/
# https://bitbucket.org/cmaurini/varfrac_for_cism
#==========================================================
#==========================================================
# Copyright (C) 2012-2013 Corrado Maurini.
# Licensed under the GNU LGPL Version 3.
#
#==========================================================
#                      Import libraries
#==========================================================
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

from dolfin import *
#import sympy
import numpy as np
from numpy import linalg as LA
from petsc4py import PETSc

try:
    from petsc4py import PETSc
except:
    print ("\n*** Warning: you need to have petsc4py installed for this module to run\n")
    print ("\nExiting.\n")
    exit()

if not has_petsc4py():
    print ("\n*** Warning: Dolfin is not compiled with petsc4py support\n")
    print ("\nExiting.\n")
    exit()

#==========================================================
#==========================================================
#                   The main function
#==========================================================
#==========================================================
def main():
    #    test_TAO_linear_bound_solver()
    #    test_SNES_linear_bound_solver()
    parameters["linear_algebra_backend"] = "PETSc"
    #==========================================================
    #                        Parameters
    #==========================================================
    
    # Some constants
    mu = 2.; p = 2.; tol = 1e-14; zero = Constant(0.0);
    
    ubVal = 0.2; lbVal = -100;
    
    # Create mesh and define function space
    Lx = 1.; Ly = 0.1; Nx = 100; Ny = 10;
    # The mesh
    #    mesh = UnitSquareMesh(Nx, Ny)
    mesh = RectangleMesh(Point(0,0), Point(Lx, Ly), Nx, Ny, "right/left") # Latest version
    
    # Degree of the varialbes
    u_degree = 2; grad_u_degree = 1; h_degree = 1;
    # Define functional spaces
    Vh = FunctionSpace(mesh, 'P', h_degree)
    Vu = FunctionSpace(mesh, 'P', u_degree)
    Vu_grad = FunctionSpace(mesh, 'P', grad_u_degree)

    # Define Dirichlet boundary conditions Rec
    def BoundaryX0(x, on_boundary):
        return on_boundary and near(x[0], 0., tol)
    
    def BoundaryX1(x, on_boundary):
        return on_boundary and near(x[0], 1., tol)
    # Apply B.C to Vu
    bcX0 = DirichletBC(Vu, zero, BoundaryX0)
    bcX1 = DirichletBC(Vu, zero, BoundaryX1)
    bcs = [bcX0, bcX1]
    
    
    #    # Define Dirichlet boundary conditions Unit Square
    #    def boundary(x, on_boundary):
    #        return on_boundary
    #    bcs = DirichletBC(Vu, zero, boundary)
    
    # Define variational problem
    u = TrialFunction(Vu)
    usol = Function(Vu)
    v = TestFunction(Vu)
    h = Function(Vh)
    
    # The RHS of the PDE
    fVal = 2.
    f = Expression('fVal', degree=1, fVal = fVal)
    
    # Initial guess for h
    h0Val = 1.
    h0 = Expression('h0Val', degree=1, h0Val = h0Val)
    h.assign(interpolate(h0, Vh))
    
    a = (h**p)*inner(grad(u), grad(v))*dx
    L = f*v*dx
    
    # Define the upper and lower bounds
    ub = interpolate(Constant(ubVal), Vu)
    lb = interpolate(Constant(lbVal), Vu)
    #-----------------------------
    # Solution with TAO (linear system formulation, working only for quadratic energies)
    #-----------------------------
    # Assemble the linear system
    A, b = assemble_system(a, L, bcs)
    # redefine the functional space for the displacement field, u
    u = Function(Vu)
    
    solver = TAOLinearBoundSolver("gpcg","gmres")
    # Set some parameters
    solver.parameters["monitor_convergence"]=True
    solver.parameters["report"]=True
    
    
    # Solve the problem
    solver.solve(A, u.vector(), b , lb.vector(), ub.vector())
    #    info(solver.parameters,True)
    
    
    # The maximum of u
    u_max = np.abs(u.vector().get_local()).max()
    print('\n u_max = \n', u_max)
    
    h_max = np.abs(h.vector().get_local()).max()
    print('\n h_max = \n', h_max)
    
    # plot the solution
    plot(u, title= "Displacement")
    plot(mesh)
    vtkfile_u = File('Output/displacement.pvd')
    vtkfile_u << u

#==========================================================
#    Find the Active set and solve the adjoint equation
#==========================================================
    u_array = u.vector().get_local()
    ub_array = ub.vector().get_local()
    # Find the active set associated to the upper bound
    ActiveSet = list(u_array[:] >(ub_array[:] - 1e-5))

    pinfty = np.inf; ninfty = -np.inf;
    
    # Define the upper and lower bounds for adjoint
    ub_adj = interpolate(Constant(pinfty), Vu)
    lb_adj = interpolate(Constant(ninfty), Vu)

    # Get upper and lower bound arrays
    ub_adj_array = ub_adj.vector().get_local();
    lb_adj_array = lb_adj.vector().get_local();
    
    #Set upper and lower bound to be zero on the ActiveSet
    for atDof, isActivePoint in enumerate(ActiveSet):
        if isActivePoint == 1 :
#            print(atDof, isActivePoint) #List the active points if needed
            ub_adj_array[atDof] = 0.; lb_adj_array[atDof] = 0.;


    #Assign values
    ub_adj.vector()[:] = ub_adj_array
    lb_adj.vector()[:] = lb_adj_array
    # Print the upper and lower arrays to see
#    print('ub_adj_array', ub_adj.vector().get_local())
#    print('lb_adj_array', lb_adj.vector().get_local())

    # Assemble the linear system
    A_adj, b_adj = assemble_system(a, L, bcs)
    
    # Define the functional space for the displacement field of the adjoint solution
    u_adj = Function(Vu)
    
    solver = TAOLinearBoundSolver("gpcg","gmres")
    # Set some parameters
    solver.parameters["monitor_convergence"]=True
    solver.parameters["report"]=True
    
    
    # Solve the problem
    solver.solve(A_adj, u_adj.vector(), b_adj , lb_adj.vector(), ub_adj.vector())
    
    # The maximum of u_adj
    u_adj_max = np.abs(u_adj.vector().get_local()).max()
    print('\n u_adj_max = \n', u_adj_max)
    
    h_max = np.abs(h.vector().get_local()).max()
    print('\n h_max = \n', h_max)

    #Print the solution of the adjoint equation
    u_adj_array = u_adj.vector().get_local()
    print('u_adj_array = ', u_adj_array)
    # plot the solution
    plot(u_adj, title= "Displacement of adj")
    plot(mesh)
    vtkfile_u_adj = File('Output/adjDisplacement.pvd')
    vtkfile_u_adj << u_adj


if __name__ == '__main__':
    main()
​
written 9 weeks ago by NhaTran108  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »