How to handle with non-linear term with power less than 1?

skystar · February 18, 2024, 9:18am

Hello,

I’m trying to solve an 2-d evolutionary non-linear PDE which has a term u^\alpha v with \alpha <1, where u,v are both variables and are non-negative in my assumption. Since it’s a non-linear term, i choose to use newton solver but it does not converge from the beginning.

Here is the information provided by using “dolfinx.log.set_log_level(log.LogLevel.INFO)”:

2024-02-18 17:02:13.172 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.173 (1923.480s) [main 2024-02-18 17:02:13.174 (1923.481s) [main 2024-02-18 17:02:13.174 (1923.482s) [main 2024-02-18 17:02:13.174 (1923.482s) [main 2024-02-18 17:02:13.174 (1923.482s) [main 2024-02-18 17:02:13.174 (1923.482s) [main 2024-02-18 17:02:13.175 (1923.482s) [main 2024-02-18 17:02:13.657 (1923.964s) [main 2024-02-18 17:02:13.659 (1923.966s) [main 2024-02-18 17:02:13.660 (1923.967s) [main 2024-02-18 17:02:13.660 (1923.967s) [main 2024-02-18 17:02:13.660 (1923.967s) [main 2024-02-18 17:02:13.660 (1923.967s) [main 2024-02-18 17:02:13.660 (1923.968s) [main 2024-02-18 17:02:13.661 (1923.968s) [main 2024-02-18 17:02:13.661 (1923.968s) [main 2024-02-18 17:02:13.661 (1923.968s) [main 2024-02-18 17:02:13.661 (1923.968s) [main 2024-02-18 17:02:13.661 (1923.968s) [main 2024-02-18 17:02:13.662 (1923.969s) [main 2024-02-18 17:02:13.662 (1923.969s) [main 2024-02-18 17:02:13.662 (1923.969s) [main ] graphbuild.cpp:395 INFO| Build local part of mesh dual graph
] ordering.cpp:203 INFO| GPS pseudo-diameter:(30) 113-14
] Topology.cpp:901 INFO| Create topology
] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 1
] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 1
] partition.cpp:233 INFO| Compute ghost indices
] MPI.cpp:99 INFO| Computing communication graph edges (using PCX algorithm). Number of input edges: 0
] MPI.cpp:156 INFO| Finished graph edge discovery using PCX algorithm. Number of discovered edges 0
] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 0
] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 0
] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 0
] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 0
] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 0
] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 0
] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 0
] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
] MPI.h:379 INFO| Number of neighbourhood source ranks in distribute_to_postoffice: 0
] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 0
] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
] MPI.h:519 INFO| Neighbourhood destination ranks from post office in distribute_data (rank, num dests, num dests/mpi_size): 0, 0, 0
] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 0
] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 0
] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 0
] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 0
] SparsityPattern.cpp:384 INFO| Column ghost size increased from 0 to 0
] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
] NewtonSolver.cpp:36 INFO| Newton iteration 2: r (abs) = inf (tol = 1e-10) r (rel) = -nan(tol = 1e-06)
] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
] NewtonSolver.cpp:36 INFO| Newton iteration 3: r (abs) = inf (tol = 1e-10) r (rel) = -nan(tol = 1e-06)
] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
] NewtonSolver.cpp:36 INFO| Newton iteration 4: r (abs) = inf (tol = 1e-10) r (rel) = -nan(tol = 1e-06)
] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
] NewtonSolver.cpp:36 INFO| Newton iteration 5: r (abs) = inf (tol = 1e-10) r (rel) = -nan(tol = 1e-06)
] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
] NewtonSolver.cpp:36 INFO| Newton iteration 6: r (abs) = inf (tol = 1e-10) r (rel) = -nan(tol = 1e-06)
] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.
] NewtonSolver.cpp:36 INFO| Newton iteration 7: r (abs) = inf (tol = 1e-10) r (rel) = -nan(tol = 1e-06)
] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.

And the following information is the same as above, r(rel) is -nan from the first step.

Then my main code:

> 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
> from dolfinx.fem import Function, functionspace
> from dolfinx.fem.petsc import NonlinearProblem
> from dolfinx.io import XDMFFile
> from dolfinx.nls.petsc import NewtonSolver
> from ufl import dx, grad, inner, dot
> 
> import dolfinx.mesh
> import dolfinx
> import gmsh
>
> dim = 2
> n_default = 8
> T_default = 1
> dt = 0.001
> num_steps_default = int(T_default / dt)
> 
> u0 = 1
> v0 = 0.1
> w0 = 1
> lamda = 1.5
> alpha = 0.5
> chi = 100
> 
> msh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, n_default, n_default, dolfinx.mesh.CellType.triangle)
> P1 = element("Lagrange", msh.basix_cell(), 1)
> ME = dolfinx.fem.functionspace(msh, mixed_element([P1, P1, P1]))
> v_1, v_2, v_3 = ufl.TestFunctions(ME)
> U = Function(ME)
> U_0 = Function(ME)
> u, v, w = ufl.split(U)
> u_last, v_last, w_last = ufl.split(U_0)

> def ic1(u_0, dim, uu, vv, ww):
>     if dim == 2:
>         u_0.sub(0).interpolate(lambda x: 1 + uu * np.exp(-1 * ((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2)))
>         u_0.sub(1).interpolate(lambda x: 1 + vv * np.exp(-1 * ((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2)))
>         u_0.sub(2).interpolate(lambda x: 1 + ww*np.exp(-1 * ((x[0]-0.5) ** 2 + (x[1]-0.5) ** 2)))
>         u_0.x.scatter_forward()

> def calculate_uh(alpha, u0, v0, w0, chi, lamda):
>
>     F1 = ((u - u_last) / dt) * v_1 * dx + dot(grad(u), grad(v_1)) * dx - chi * dot(u * grad(v), grad(
>         v_1)) * dx - lamda * v_1 * dx + u * v_1 * dx + (u ** alpha) * w * v_1 * dx
>     F2 = ((v - v_last) / dt) * v_2 * dx + dot(grad(v), grad(v_2)) * dx + v * v_2 * dx - (u ** alpha) * w * v_2 * dx
>     F3 = ((w - w_last) / dt) * v_3 * dx + dot(grad(w), grad(v_3)) * dx + w * v_3 * dx - v * v_3 * dx
>     F = F1 + F2 + F3
> 
>     problem = NonlinearProblem(F, U)
>     solver = NewtonSolver(MPI.COMM_WORLD, problem)
>     solver.convergence_criterion = "incremental"
>     solver.rtol = 1e-6
>     solver.report = True
> 
>     ksp = solver.krylov_solver
>     opts = PETSc.Options()
>     option_prefix = ksp.getOptionsPrefix()
>     opts[f"{option_prefix}ksp_type"] = "preonly"
>     opts[f"{option_prefix}pc_type"] = "lu"
>     ksp.setFromOptions()
> 
>     ic1(U_0, dim, u0, v0, w0)
>     t = 0
>     for i in range(num_steps_default):
>         t = t + dt
>         try:
>             solver.solve(U)
>             U_0.x.array[:] = U.x.array
>         except:
>             print('Blow up. t = '+str(t))
>             return U_0
>     return U
>
> dolfinx.log.set_log_level(log.LogLevel.INFO)
> U = calculate_uh(alpha, u0, v0, w0, chi, lamda)

Besides, the situations when \alpha \geq 1 and \alpha =0 worked, so i think the problem is about newton solver and the non-linear term.

Thanks so much for providing any advice!

francesco-ballarin · February 18, 2024, 11:10am

Please make sure to include the entire code in a block created with “```”, and which does not contain “>”, otherwise it’s difficult to copy, paste and run your code.

From what I get by reading your code, you have a time loop, and you assign the initial condition to U_0. However, I can’t see where you also assign U_0 to U, i.e. where you say that at time dt the nonlinear solver should start from a solution U which has the same values as U_0 as initial guess. If indeed you do not do that, the nonlinear solver will start from the zero solution.

If that is indeed the issue, note that it will only affect the very first time step. For subsequent one, the solver will use whatever value it finds in U as initial guess. In particular, the actual value it finds will be the converged solution at the previous time stemp, unless you do something to change the content of U.

skystar · February 18, 2024, 11:58am

Sorry for the code, but your advice really works, thanks!

Topic		Replies	Views
Non-linear solver convergence failure dolfinx	1	535	December 14, 2021
Non linear solver fails to converge. Is there a thumb-rule to avoid this in the future? variational formulation	7	2627	November 4, 2021
No convergence for nonlinear elastic problem with spectral decomposition constitutive equation dolfinx	1	338	April 8, 2024
1D nonlinear solver does not converge	0	484	August 25, 2020
Is my problem linear? Newton solver converges in 5 steps mathematics dolfinx , newton	2	258	December 16, 2023

How to handle with non-linear term with power less than 1?

Related topics