179
views
0
4 months ago by
Greetings!
I have an iterative process:
1. Get state from nonlinear boundary problem.
2. Solve conjugate system to get gradient.
3. Recalculate control function.

Problem is quite simple: time for the first step grows dramatically with every iteration:
   1: 0.0585370063782
2: 0.0657551288605
3: 0.0835440158844
4: 0.0911159515381
....
99: 1.99743199348
100: 1.69048595428
101: 1.64172911644
102: 1.69205212593
103: 1.65970301628​

In these conditions calculate, for example, a thousand of iterations might be impossible. Need to mention that a number of newton iterations remains the same (~2--3). So I have no any idea why that's happening.

Code to reproduce:

from __future__ import print_function
import time as tm
from dolfin import *
from dolfin.cpp.common import set_log_level
from dolfin.cpp.mesh import UnitCubeMesh

set_log_level(50)
def main():
omega = UnitCubeMesh(5,5,5)
a, alpha, b, epsilon = 0.06, 0.3333, 0.025, 0.000000000001
finite_element = FiniteElement("CG", omega.ufl_cell(), 1)
state_space = FunctionSpace(omega, finite_element * finite_element)
state = Function(state_space)
theta, phi = split(state)
v, h = TestFunctions(state_space)
theta_0 = interpolate(Constant(1), FunctionSpace(omega, finite_element))
u_control = Expression('0.3', degree=3)
theta_b = 0

def set_and_solve_boundary():
+ b * inner(theta ** 4 - phi, v) * dx \
+ inner(phi - theta ** 4, h) * dx \
+ inner(theta - theta_b, v) * ds - u_control * h * ds
solve(Boundary_problem == 0, state)
return state.split()

def set_and_solve_conjugate():
p_1, p_2 = TrialFunctions(state_space)
+ 4 * theta ** 3 * inner(b * p_1 - p_2, v) * dx \
+ inner(p_1, v) * ds \
+ inner(p_2 - b * p_1, h) * dx
J_theta = - (theta - theta_0) * v * ds

for _ in range(1000):
start = tm.time()
set_and_solve_boundary()
print(str(_).rjust(4, " ")+": "+str(tm.time() - start))
p_1, p_2 = set_and_solve_conjugate()
u_control = Expression('u + (p_2 - eps*u)', u=u_control,
p_2=p_2, eps=epsilon, degree=3)
main()
​

I will appreciate any suggestions
Community: FEniCS Project

 u_control = interpolate(Expression('u + (p_2 - eps*u)', u=u_control,
​