Time dependent non linear SNES solver

dolfinx
1

Hello,

Are there any changes necessary for to the code provided in

to solve a time-dependent non-linear problem?

Perhaps with these following lines?

b = la.create_petsc_vector(V.dofmap.index_map, V.dofmap.index_map_bs)
J = create_matrix(problem.a)

I tried with creating a NonlinearPDE_SNESProblem instance inside a loop, since F for my problem depends on the previous step solution. But the solution doesn’t change with time. This makes me suspect that something needs updating while running in a loop.

Can someone please point me in the right direction?

Thank you
Krishna Kenja

2

It would be useful of you could supply a minimal example illustrating the issue you are having

3

Hello @dokken

I have a problem, I am trying to solve a nonlinear problem. What happens is that after a few moments of time, it gives a message that it does 0 iterations and the solution does not change.

import fenics as fe
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
def model2(t, a, b, c, d, e, f):
    return a*np.exp(-((t-b)/c)**2) + d*np.exp(-((t-e)/f)**2)
TIME_STEP_LENGTH = 0.05
SIMULATION_TIME = 3600*23
N_TIME_STEPS = int(SIMULATION_TIME/TIME_STEP_LENGTH)
t = np.linspace(0, SIMULATION_TIME, N_TIME_STEPS)
Nx = Ny = 20
mesh = UnitSquareMesh(Nx, Ny)
class Inflow(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0)
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0)
inflow = Inflow()
bot = Bottom()
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
w = fe.TestFunction(Q)
T_t0 = fe.Constant(296.35)
u0 = fe.Expression(('0.25', '0.0'), degree=2)
u = fe.interpolate(u0, V)
T_prev = fe.interpolate(T_t0, Q)
T_next = fe.Function(Q)
gamma = 0.01
sigma_B = 5.6703 * 1e-8
T0 = fe.Constant(292.15)
T0_prev = fe.interpolate(T0, Q)
T_inflow = fe.Expression('T0 - gamma*x[1]', degree=1, T0=T0, gamma=gamma)
bc_in = DirichletBC(Q, T_inflow, inflow)
bc_bot = DirichletBC(Q, fe.Constant(292.15), bot)
bc = [bc_in, bc_bot]
rho = 1.1614    # density (kg m-3)
mu = 1.87 * 1e-5    # dynamic viscosity (Pa s)
cp = 1005   # specific heat (J kg-1 K-1)
k0 = fe.Constant(0.375)
e0 = fe.Constant(0.075)
k_prev = fe.interpolate(k0, Q)
e_prev = fe.interpolate(e0, Q)
mut = rho * 0.09 * k_prev ** 2 / e_prev
t_data = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]) * 3600
ft_data = np.array([1.60, 1.60, 1.50, 1.60, 1.50, 1.40, 1.80, 46.10, 199.30, 375.00, 485.30, 576.80, 538.80, 514.30,
                    316.80, 262.30, 136.50, 92.70, 10.40, 1.40, 1.40, 1.30, 1.50, 1.40])
initial_point = [200, 10, 2, 500, 15, 3]
x, _ = curve_fit(model2, t_data, ft_data, p0=initial_point)
Rs = model2(t, *x)
n=0
Q_rad = (1 - 0.27) * Rs[n] + 0.77 * sigma_B * T_next ** 4 - 0.85 * sigma_B * T0_prev ** 4
F = (rho * cp * 1.0 / TIME_STEP_LENGTH * fe.dot((T_next - T_prev), w) * fe.dx
     + rho * cp * fe.dot(fe.dot(u, fe.grad(T_next)), w) * fe.dx
     + (mu + mut) * fe.inner(fe.grad(T_next), fe.grad(w)) * fe.dx
     - 1.0 / 0.1229 * fe.dot(Q_rad, w) * fe.dx)
for n in range(N_TIME_STEPS):
    solve(F == 0, T_next, bc, solver_parameters={"newton_solver": {"relative_tolerance": 1e-6}})
    T_prev.assign(T_next)
4

I guess that means that your problem reaches a steady-state?
You could also experiment with different linear solvers ({"linear_solver":"mumps"} inside solver_parameters for instance), and visualize the output to see if you are getting the expected solution.