### Non-linear coupled equations are not converging

197

views

1

Hello everyone

I am trying to simulate non-linear coupled equations and this is my code:

The problem is with changing coefficients or time, either is not converging or it is converging to non-sense(value of u_1 for example becomes zero in whole domain).

What is problem in your opinion?!

Initial condition seems reasonable. Sometimes I think problem is variational formulation, but I am not sure.

These are the equations:

$\frac{\partial u_1}{\partial t}=D\times\left(\Delta u_1+\frac{K}{u_2}\times\Delta u_2\right)\backslash\frac{\partial u_2}{\partial t}-\left(\frac{K}{c}\right)\times\left(\frac{u_2}{u_1}\right)\times\frac{\partial u_1}{\partial t}=X\times\Delta u_2$∂

I am trying to simulate non-linear coupled equations and this is my code:

```
# -*- coding: utf8 -*-
from __future__ import print_function
from fenics import *
import numpy as np
from dolfin import *
from mshr import *
tol = 1E-14
t = 0.0
class InitialConditions(Expression):
def eval(self, values, x):
values[0] = 0.0
values[1] = 315
if (near(x[0], 2.0, 0.2) or near(x[0], 13.0, 0.2)) and \
near(x[1], 1.0, 0.1):
values[0] = 5000.0
values[1] = 315.01
def value_shape(self):
return (2,)
T = 1.0e-1 # final time
num_steps = 100 # number of time steps
dt = T / num_steps # time step size
D = 1.000e-10 # Diffusion Coefficient
K = 8.000e-1 # Thermometric diffusion ratio
c = 4.000e3 # Heat Capacity
X = 2.000e-7 # Thermoemtric diffusion coefficient
# Create mesh and define function space
mesh = RectangleMesh(Point(0., 0.),Point(15.0, 2.0), 30, 30)
P1 = FiniteElement("P", triangle,2)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)
# Define test functions
v_1, v_2 = TestFunctions(V)
# Define functions for concentration and Temperature
u = Function(V)
u_n = Function(V)
# Split system functions to access components
u_1, u_2 = split(u)
# Define boundary condition
z_1 = Expression(('0.0','315.0'), degree = 2)
z_2 = Expression(('0.0','315.0'), degree = 2)
#tol = 1E-14
def boundary_D(x, on_boundary):
return on_boundary and not ((near(x[0], 2.0, 0.5) or near(x[0], 13.0, 0.5)) and\
near(x[1], 2.0, tol)) and not (near(x[0], 0.0, tol) or\
near(x[0], 15.0, tol))
def boundary_D_source(x, on_boundary):
return on_boundary and (near(x[0], 2.0, 0.5) or near(x[0], 13.0, 0.5)) and \
near(x[1], 2.0, tol)
bc_1= DirichletBC(V, z_1, boundary_D_source)
bc_2 = DirichletBC(V, z_2, boundary_D)
bcs = [bc_1, bc_2]
#bcs = bc_2
# Define initial value
f = InitialConditions(degree = 1)
u_n = interpolate(f, V)
S_1 = Source(degree = 2)
u_n1, u_n2 = u_n.split()
# Define variational problem
v_1, v_2= TestFunctions(V)
k = Constant(dt)
K = Constant(K)
D = Constant(D)
c = Constant(c)
X = Constant(X)
g = Constant(-1.0)
r = Constant(-1.00e-6)
F = u_2*((u_1 - u_n1))*v_1*dx + k*D*dot(grad(u_1), grad(v_1*u_2))*dx \
+ k*D*g*v_1*u_2*ds + k*D*K*dot(grad(u_2), grad(v_1))*dx + k*D*K*r*v_1*ds\
+ u_1*((u_2 - u_n2))*v_2*dx - (K/c)*(u_2)*((u_1 - u_n1))*v_2*dx\
+ k*X*dot(grad(u_2), grad(v_2*u_1))*dx + k*X*r*v_2*u_1*ds
vtkfile_1 = File("Project/u1.pvd")
vtkfile_2 = File("Project/u2.pvd")
du = TrialFunction(V)
# Time-stepping
for n in range(num_steps):
_u_1, _u_2 = u_n.split()
vtkfile_2 << (_u_1,t)
# Compute solution
J = derivative(F, u, du)
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()
# Update previous solution
u_n.assign(u)
t += dt
```

The problem is with changing coefficients or time, either is not converging or it is converging to non-sense(value of u_1 for example becomes zero in whole domain).

What is problem in your opinion?!

Initial condition seems reasonable. Sometimes I think problem is variational formulation, but I am not sure.

These are the equations:

$\frac{\partial u_1}{\partial t}=D\times\left(\Delta u_1+\frac{K}{u_2}\times\Delta u_2\right)\backslash\frac{\partial u_2}{\partial t}-\left(\frac{K}{c}\right)\times\left(\frac{u_2}{u_1}\right)\times\frac{\partial u_1}{\partial t}=X\times\Delta u_2$∂

`u`_{1}∂`t`=`D`×(Δ`u`_{1}+`K``u`_{2}×Δ`u`_{2})\∂`u`_{2}∂`t`−(`K``c`)×(`u`_{2}`u`_{1})×∂`u`_{1}∂`t`=`X`×Δ`u`_{2}
Community: FEniCS Project

I cannot reproduce your problem since I get a different error in line 70 NameError: name 'Source' is not defined.

written
5 months ago by
Eleonora Piersanti

Oh I am sorry. Delete that line. I forgot to delete it, it does not affect the code. Thank you for checking my code by the way.

written
5 months ago by
Morteza

### 1 Answer

0

When dealing with such low values for coefficients such as the diffusion and thermometric diffusion maybe it would be a good idea to nondimensionalize your equation. We need to be careful with too little and too big values as we are dealing with float points arithmetic.

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