### How to implement different initial conditions on different subdomains

71
views
0
19 days ago by
Hi
I have a diffusion problem. I used to have only 1 variable (c1) and a single domain. The variable has an initial value (X1) on the whole domain:

That is a part of the code and how I used to solve it:

from fenics import *
import numpy as np

#Defining Mesh
mesh = RectangleMesh(Point(0,0), Point(1, 1), 4, 4, "right")

#Defining Element
Element1 = FiniteElement("CG", mesh.ufl_cell(), 1)

# Defining the mixed function space
W_elem = MixedElement([Element1])

W = FunctionSpace(mesh, W_elem)

z = Function(W)
dz=TrialFunction(W)

#Defining the test functions
v_1 = TestFunctions(W)

#Boundary conditions
bcs = [BCs]

# Defining  function spaces for the variable c1
V_c = FunctionSpace(mesh, Element1)

# Initial Condition
C_INITIAL = interpolate(X1, V_c)

#Part of the variational form
F = c1 * v_1 * dx + dt * dot(grad(c1), grad(v_1))*dx

#Solve
while t <= T:
J = derivative(F, z, dz)
problem = NonlinearVariationalProblem(F, z, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()

C_INITIAL.assign(c1)  #Assigning INITIAL condition

t += dt
f << c1

And everything would go well. Now I want to solve a different problem where I have 3 variables (c1 , c2 and c3)  instead of 1. In addition my domain has been splitted into 2 subdomains and initial values of the variables are different in these two subdomains:

The modified code is:

from fenics import *

#Defining Mesh
mesh = RectangleMesh(Point(0,0), Point(1, 1), 4, 4, "right")

class A (SubDomain):
def inside(self, x, on_boundary):
return True if x[0] >= 0.0 and x[0] <= 0.5 and (x[1] >= 0.0 ) and (x[1] <= 1.0 ) else False

class B (SubDomain):
def inside(self, x, on_boundary):
return True if x[0] >= 0.5 and x[0] <= 1.0 and (x[1] >= 0.0 ) and (x[1] <= 1.0 ) else False

domains = CellFunction("size_t", mesh)
domains.set_all(0)

left = A ()
left.mark(domains, 1)

right = B()
right.mark(domains, 2)

#Defining Element
Element1 = FiniteElement("CG", mesh.ufl_cell(), 1)

# Defining the mixed function space
W_elem = MixedElement([Element1,Element1,Element1 ])

W = FunctionSpace(mesh, W_elem)

z = Function(W)
dz=TrialFunction(W)

#Splitting the variables
c1 , c2, c3 = split(z)

#Defining the test functions
v_1,v_2,v_3  = TestFunctions(W)

#Boundary conditions
bcs = [BCs]

# Defining  function spaces for the variable c1
V_c = FunctionSpace(mesh, Element1)

# Initial Condition
#How to define and interpolate different initial values for 3  variables in subdomains?

#Part of the variational form
F1 = c1 * v_1 * dx(1) + dt * dot(grad(c1), grad(v_1))* dx(2)
F2 = c2 * v_2 * dx(1) + dt * dot(grad(c2), grad(v_2))* dx(2)
F3 = c3 * v_3 * dx(1)  + dt * dot(grad(c3), grad(v_3))* dx(2)

F = F1 + F2 + F3

#Solve
while t <= T:
J = derivative(F, z, dz)
problem = NonlinearVariationalProblem(F, z, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()

(c1, c2, c3) = z.split(True)

#How to assign different initial conditions in different domains for 3 variable??

t += dt
f << c1

Does anybody know how I can implement different initial conditions for 3 different variables in 2 different subdomains?
Thanks!
Community: FEniCS Project

2
17 days ago by

I notice that there is an overlap of the domains A and B at x = 0.5. Assuming that the division of the subdomains is a vertical line at x=0.5 (and assuming 0.5 belongs to A), i would do something like this:

class initExpression(Expression):
def __init__(self, **kwargs):
self.boundary = kwargs["boundary"]
self.x1 = kwargs["x1"]
self.x2 = kwargs["x2"]
self.y1 = kwargs["y1"]
self.y2 = kwargs["y2"]
self.z1 = kwargs["z1"]
self.z2 = kwargs["z2"]

def eval(self, value, x):
value[0] = self.x1
value[1] = self.y1
value[2] = self.z1
if x[0] > self.boundary:
value[0] = self.x2
value[1] = self.y2
value[2] = self.z2

and to initialize it:
assuming boundary at x = 0.5, x1 = y1 = z1 = 1  and x2 = y2 = z2 = 2

u_init = initExpression(boundary = 0.5, x1 = 1.0, x2 = 2.0, y1 = 1.0, y2 = 2.0, z1 = 1.0, z2 = 2.0, element=W.ufl_element())
z.assign(u_init)
Thanks for your reply. Although it did not work. As you can see we should assign the initial condition to the variable that we are calculating (This variable is c in my problem) in the solver but before that, we need to interpolate the initial condition on the FunctionSpace. If we use :
C_previous = interpolate(u_init, V_c)
​
This error pops up:

*** Error:   Unable to interpolate function into function space.
*** Reason:  Rank of function (1) does not match rank of function space (0).​
Based on what you suggested, I simplified the problem to reduce the number of the variables to only 1 (c).  I took the mid boundary at x = 0.005 and the value of the initial conditions constant on both domain. In other words I am using only 1 initial condition on the whole domain but I am just splitting the domain to 2 different domains and applying the initial condition separately on these domain. It is obvious that the same results is expected to obtain as there was only 1 domain and 1 initial condition on the whole domain. Here is what I did based on your suggestion:

class initExpression(Expression):
def __init__(self, **kwargs):
self.boundary = kwargs["boundary"]
self.x1 = kwargs["x1"]
self.x2 = kwargs["x2"]

def eval(self, value, x):
value[0] = self.x1

if x[0] > self.boundary:
value[0] = self.x1

else:
value[0] = self.x2

C_0 = Constant(500.0)

u_init = initExpression(boundary = 0.005, x1 = C_0, x2 = C_0, element=W.ufl_element())

C_previous = interpolate(u_init, V_c)
​

This is what I did in the solver:

while t <= T:
J = derivative(F, z, dz)
problem = NonlinearVariationalProblem(F, z, bcs, J)
solver = NonlinearVariationalSolver(problem)

solver.solve()

# Assign the initial condition to the variable "c"
C_previous.assign(c)

t += dt​
And the same error pops up. Any better idea how to implement it? Thanks
written 14 days ago by Leo
1
With a little modification, I resolved this issue:

class initExpression(Expression):
def __init__(self, **kwargs):
self.boundary = kwargs["boundary"]
self.x1 = kwargs["x1"]
self.x2 = kwargs["x2"]

def eval(self, value, x):

if x[0] > self.boundary:
value[0] = self.x1
else:
value[0] = self.x2

u_init = initExpression(boundary = 0.005, x1 = C_0, x2 = C_0, degree=0)
C_previous = interpolate(u_init, V_c)
.
.
.
#Solve

C_previous.assign(c)​
written 14 days ago by Leo