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

71

views

0

Hi

I have a diffusion problem. I used to have only

Does anybody know how I can implement different initial conditions for 3 different variables in 2 different subdomains?

Thanks!

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

### 1 Answer

2

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)
```

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

Please login to add an answer/comment or follow this question.

This error pops up:

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:This is what I did in the solver:

And the same error pops up. Any better idea how to implement it? Thanks