Implementation of boundary conditions from file

Asked by Bjarnthor Egilsson on 2013-10-25

Hi, I'm doing a simulation where I want to read a velocity field from a stored solution and implement that solution as a boundary condition on a different mesh., using CBC Fully segregated solver.

Here is what I have so far:
    def create_boundaries(self):
        OldMesh = Mesh('geometry/OldMesh.xml')
        V_old = VectorFunctionSpace(OldMesh,'CG',1)
        uxml = 'results/xml/u_1.xml.gz'
        uold = Function(V_old,uxml)
        V = VectorFunctionSpace(self.mesh,'CG',1)
        unew = project(uold,V)

Now I know this boundary condition can be specified by:
DirichletBC(V,unew,boundary_ID), but my problem is how can I apply DirichletBCs directly using cbc? And do I have to do something special for the segregated solver, like decomposing the velocity into components?

Structure of the code (stripped down for clarification):

from cbc.pdesys import *
from NSProblem import *

class VR(NSProblem):

    def __init__(self, parameters):
        NSProblem.__init__(self, parameters=parameters)
        self.mesh = Mesh('geometry/Mesh.xml')
        self.boundaries = self.create_boundaries()
        hmin = MPI.min(self.mesh.hmin())
        problem_parameters['dt'] = cfl * hmin**2 /(Umax*(nu+hmin*Umax))
        # To initialize solution set the dictionary q0:
        self.q0 = initial_velocity
    def create_boundaries(self):
        OldMesh = Mesh('geometry/Combined.xml')
        V_old = VectorFunctionSpace(OldMesh,'CG',1)
        uxml = 'results/xml/u_1.xml.gz'
        uold = Function(V_old,uxml)
        V = VectorFunctionSpace(self.mesh,'CG',1)
        unew = project(uold,V)
        plot(unew)
        bnd = MeshSubDomain(3, 'Wall') #Here should the boundary condition be specified

        return [bnd]

    def prepare(self):
        #Update the bc here, at each time step

if __name__ == '__main__':
    from cbc.cfd.icns import NSFullySegregated, NSSegregated, solver_parameters
    import time
    parameters["linear_algebra_backend"] = "PETSc"
    problem_parameters['viscosity'] = nu
    problem_parameters['T'] = 2
    problem_parameters["time_integration"] = 'Transient'

    solver_parameters = recursive_update(solver_parameters,
    dict(degree=dict(u=2,u0=1,u1=1),
         pdesubsystem=dict(u=101, p=101, velocity_update=101),
         linear_solver=dict(u='gmres', p='gmres', velocity_update='gmres'),
         precond=dict(u='jacobi', p='hypre_amg', velocity_update='jacobi'))
         )
    problem = VR(problem_parameters)
    solver = NSFullySegregated(problem, solver_parameters)
    problem.solve()

Question information

Language:
English Edit question
Status:
Solved
For:
CBC.PDESys Edit question
Assignee:
No assignee Edit question
Solved by:
Mikael Mortensen
Solved:
Last query:
Last reply:

Hi,

You need to create an expression, attach the function to be called at the boundary and overload eval. If you look at the turbulence demos, all turbulence problems use the solution from a turbulent channel as inlet condition. Basically:

class My_inlet_ux(Expression)
    def __init__(self, u_old):
        self.u_old = u_old

    def eval(self, values, x):
        values[0] = self.u_old(x)[0]

my_inlet_ux = My_inlet(u_old)

and then similarily for the other components. Hope you can work out the details:-)

Mikael

Thanks Mikael Mortensen, that solved my question.