Hello everyone, I am interested in solving a tutorial problem in fenicsx with time increment adaptive based on iteration changes or number of times solver is solving. So as a general logic, I am writing the same in the code like,
while (t < T):
t += float(dt)
r = solver.solve(u)
file.write_function(c, t)
ii += 1
print(f"Step {int(ii)}: num iterations: {r[0]}")
u.x.scatter_forward
u0.x.array[:] = u.x.array
if 10 < ii <= 100 :
dt = 1e-8
elif 100 < ii <= 1000:
dt = 1e-7
elif 1000 < ii <= 10000:
dt = 1e-6
Like above when I implemented, I observe the time increment happening and the solver is converging every time but unfortunately the output related to running with constant time step yields some different result. So could anyone please confirm me whether these sort of application of adaptive time increment is possible or this method of implementation is wrong?. Here it is implemented based on number of times solving, but in addition I tried based on number of iterations as well. All these not yielding exact result as doing the same with constant timestep. Is there any specific ways to do some time increment changes other than direct implementation like above?
Any advice on the above is greatly appreciated. Thanks everyone in advance.
Adaptive time stepping is entirely possible. However, in your proposed implementation it’s unclear how dt is being used. Typically you’d wrap dt in a Constant such that its value may be modified prior to reassembly of your system for computation of the next time step. Some pseudocode would be:
...
dt = dolfinx.fem.Constant(mesh, 1.0)
...
a = ufl.inner((u - un) / dt, v) * ufl.dx + ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
...
while t < T:
...
dt.value = some_new_dt_value_based_on_logic
...
dolfinx.fem.petsc.assemble_matrix(A, a, bcs=bcs)
A.assemble()
...
# Solve next time step
...
Thanks for your response sir. I am exactly using it in the same as in Cahn Hilliard tutorial, but when I noticed in my whole code,
import os
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, plot,mesh
from dolfinx.fem import Function, functionspace,assemble_scalar
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_unit_square
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner
import dolfinx
lmbda = 1.0e-02 # surface parameter
msh = create_unit_square(MPI.COMM_WORLD, 1, 1, CellType.triangle)
P1 = element("Lagrange", msh.basix_cell(), 1)
ME = functionspace(msh, mixed_element([P1, P1]))
q, v = ufl.TestFunctions(ME)
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
# Split mixed functions
c, mu = ufl.split(u)
c0, mu0 = ufl.split(u0)
##################################
dt = 1e-7
dt1 = dolfinx.fem.Constant(msh, dt) # time step
################################## Modification as per suggestion
# Zero u
u.x.array[:] = 0.0
# Interpolate initial condition
u.sub(0).interpolate(lambda x: 0.63 + 0.01 * (0.5 - np.random.rand(x.shape[1])))
u.x.scatter_forward()
# Compute the chemical potential df/dc
c = ufl.variable(c)
f = 25 * c**2 * (1 - c)**2
dfdc = ufl.diff(f, c)
# Weak statement of the equations
F0 = inner(c, q) * dx - inner(c0, q) * dx + dt * inner(grad(mu), grad(q)) * dx
F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
F = F0 + F1
# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-5
# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options() # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}absolute_tolerance"] = 1e-8
opts[f"{option_prefix}relative_tolerance"] = 1e-5
ksp.setFromOptions()
# Step in time
t = 0.0
# Reduce run time if on test (CI) server
if "CI" in os.environ.keys() or "GITHUB_ACTIONS" in os.environ.keys():
T = 3 * dt
else:
T = 1000* dt
# Get the sub-space for c and the corresponding dofs in the mixed space
# vector
V0, dofs = ME.sub(0).collapse()
c = u.sub(0)
u0.x.array[:] = u.x.array
ii=0
while (t < T):
ii += 1
t += float(dt1)
r = solver.solve(u)
print(f"Step {int(t/dt)}: num iterations: {r[0]}")
u0.x.array[:] = u.x.array
if ii > 500:
dt = 1e-5
By adding like the above can modify my dt in whole weak form while solving or am I understood anything wrongly? Could you please confirm me once. I am doubting that changing dt will create impact while solving weak form in solver or just increasing in time. As previously showed code not created impact in output. So i was confused. Thanks again for your help sir.
Thanks a lot sir, it works. Currently the results are different and seems to be case where solver running with changed timestep. Additionally one more question, sometimes my solver is exceeded with newton iterations and stops the code. I am aiming to reduce timestep if my solver is taking more than say 25 iterations in Newton solver. With the above code, may i able to do that with addition in while loop like,
else:
print ("Newton iteration not converged, reduce dt")
bisection_count += 0
if bisection_count > 5:
print("Error: Too many bisections")
break
print("Error Halfing Time Step")
t = t - dt
dt = dt / 2
print(f"New Time Step: {dt}")
u.x.array[:] = u0.x.array
I am trying to do these, but since my newton solver directly coming out if newton solver exceeding the number of iterations when 50 reached and stops the code. So may I know is it possible to make above change only with Custom newton solver or I can able to add these in above code where the NewtonSolver is used by importing. Can you please some help over the above query sir. Thanks for your continuous support.