Different result in FeniCSx from the 'almost same' legacy FeniCS code (problem related to dolfinx_mpc?)

dolfinx
, ,
1

I’m using the latest dolfinx version 0.10.0.0, and I’m trying to replicate the FeniCS code (Implementing k-epsilon model in FEniCS - how to implement terms in variational formulation that depend on a solution?) but in FeniCSx (a fully developed channel flow using a Lam-Bremhorst k-ε model with periodic BC in x-direction.)

I think that the problem might be related to how I’m implementing the loop, and to the periodic BC that I’m imposing with dolfinx_mpc.

This is the part of the code where I’m defining the mpc constraints:

# Create multi point constraint
mpc_u = MultiPointConstraint(V)
mpc_u.create_periodic_constraint_topological(V, ft, 2, periodic_relation, bcu)
mpc_u.finalize()

mpc_p = MultiPointConstraint(Q)
mpc_p.finalize()

mpc_k = MultiPointConstraint(K)
mpc_k.create_periodic_constraint_topological(K, ft, 2, periodic_relation, bck)
mpc_k.finalize()

mpc_e = MultiPointConstraint(K)
mpc_e.create_periodic_constraint_topological(K, ft, 2, periodic_relation, bce)
mpc_e.finalize()


# Define trial and test functions
u = TrialFunction(V)
p = TrialFunction(Q)
k = TrialFunction(K)
e = TrialFunction(K)

v = TestFunction(V)
q = TestFunction(Q)
phi = TestFunction(K)
psi = TestFunction(K)

# Define functions for quantities at the previous time-step:
# (also giving an inital value)
u0 = Function(mpc_u.function_space)
u0_expr = lambda x: np.vstack((u_init * np.ones_like(x[0]), np.zeros_like(x[0])))
u0.interpolate(u0_expr)
u0.x.scatter_forward()
k0 = Function(mpc_k.function_space)
k0.interpolate(lambda x: k_init * np.ones_like(x[0]))
k0.x.scatter_forward()
e0 = Function(mpc_e.function_space)
e0.interpolate(lambda x: e_init * np.ones_like(x[0]))
e0.x.scatter_forward()
# Define functions where solutions will be stored (for the new time-step):
u1 = Function(mpc_u.function_space)
p1 = Function(mpc_p.function_space)
k1 = Function(mpc_k.function_space)
e1 = Function(mpc_e.function_space)

While, this is the loop that is used to solve the equations:


# main loop
iter_max = 2 # 1000 (or: 2000# Remembering that: dt = 0.05 (1000 * 0.05 sec = 50 sec <--- Tempo totale simulato)
tol = 1e-5
print_max = False  
for iter in range(iter_max):
    # N-S in 3 steps (Chorin method):
    F1  = (1. / dt) * dot(u - u0, v)*dx \
        + dot(dot(u0, nabla_grad(u)), v)*dx \
        + (nu + nut) * inner(grad(u), grad(v))*dx 

    a_1 = form(lhs(F1)); l_1 = form(rhs(F1))
    problem1 = dolfinx_mpc.LinearProblem(a_1, l_1, mpc_u, bcs=bcu, petsc_options=petsc_options)
    u1 = problem1.solve()
    assert isinstance(u1, Function)
    if print_max:
        u1_values = u1.x.array[:]
        print(f'np.max(u1): {np.max(u1_values)}')
        print(f'np.min(u1): {np.min(u1_values)} \n')

    a_2 = dot(grad(p), grad(q))*dx; l_2 = (-1. / dt) * dot(div(u1), q)*dx 
    problem2 = dolfinx_mpc.LinearProblem(a_2, l_2, mpc_p, bcs=bcp, petsc_options=petsc_options)
    p1 = problem2.solve()
    assert isinstance(p1, Function)
    if print_max:
        p1_values = p1.x.array[:]
        print(f'np.max(p1): {np.max(p1_values)}')
        print(f'np.min(p1): {np.min(p1_values)} \n')

    a_3 = form(dot(u, v)*dx); l_3 = form(dot(u1, v)*dx - dt * dot(grad(p1), v)*dx)
    problem3 = dolfinx_mpc.LinearProblem(a_3, l_3, mpc_u, bcs=bcu, petsc_options=petsc_options)
    u1 = problem3.solve()
    assert isinstance(u1, Function)
    if print_max:
        u1_values = u1.x.array[:]
        print(f'np.max(u1): {np.max(u1_values)}')
        print(f'np.min(u1): {np.min(u1_values)} \n')

    # k-epsilon tubolence modeling:
    # Define helping functions for k-eps model
    Re_t = (1. / nu) * (k0**2 / e0)
    Re_k = (1. / nu) * (sqrt(k0) * dfun)

    f_nu = Max(Min((1 - ufl.exp(- 0.0165 * Re_k))**2 * (1 + 20.5 / Re_t), \
                Constant(mesh, 1.0)), Constant(mesh, 0.01116))
    f_1  = 1 + (0.05 / f_nu)**3 
    f_2  = 1 - ufl.exp(- Re_t**2) 
    S_sq = 2 * inner(ufl.sym(nabla_grad(u1)), ufl.sym(nabla_grad(u1)))

    # Define terms that appear in variational formulation for k-eps model
    nut     = 0.09 * f_nu * (k0**2 / e0)
    prod_k  = nut * S_sq
    react_k = e0 / k0
    prod_e  = 1.44 * react_k * f_1 * prod_k
    react_e = 1.92 * react_k * f_2

    FK  = (1. / dt) * dot(k - k0, phi)*dx \
        + dot(dot(u1, nabla_grad(k)), phi)*dx \
        + (nu + nut / 1.0) * inner(grad(k), grad(phi))*dx \
        - dot(prod_k, phi)*dx + dot(react_k * k, phi)*dx 

    a_k = form(lhs(FK)); l_k = form(rhs(FK))
    problem_k = dolfinx_mpc.LinearProblem(a_k, l_k, mpc_k, bcs=bck, petsc_options=petsc_options)
    k1 = problem_k.solve()
    assert isinstance(k1, Function)
    # print(f'k-equation --> solved!')
    if print_max:
        k1_values = k1.x.array[:]
        print(f'np.max(k1): {np.max(k1_values)}')
        print(f'np.min(k1): {np.min(k1_values)} \n')


    FE  = (1. / dt) * dot(e - e0, psi)*dx \
        + dot(dot(u1, nabla_grad(e)), psi)*dx \
        + (nu + nut / 1.3) * inner(grad(e), grad(psi))*dx \
        - dot(prod_e, psi)*dx + dot(react_e * e, psi)*dx 

    a_e = form(lhs(FE)); l_e = form(rhs(FE))
    problem_e = dolfinx_mpc.LinearProblem(a_e, l_e, mpc_e, bcs=bce, petsc_options=petsc_options)
    e1 = problem_e.solve()
    assert isinstance(e1, Function)
    # print(f'epsilon-equation --> solved!')
    if print_max:
        e1_values = e1.x.array[:]
        print(f'np.max(e1): {np.max(e1_values)}')
        print(f'np.min(e1): {np.min(e1_values)} \n')
        print(f'End of iteration: {iter}')


    # Make sure k and eps are positive
    low_bound = 1e-16
    dimension = len(k1.x.array)
    k1.x.array[:] = np.max([k1.x.array[:], low_bound * np.ones(dimension)], axis=0)  
    e1.x.array[:] = np.max([e1.x.array[:], low_bound * np.ones(dimension)], axis=0)     

    # Check if converged
    err = u0 - u1
    error_form = form(dot(err, err) * dx)
    error = np.sqrt(assemble_scalar(error_form))
    if error < tol:
        print('Fully developed flow has been reached in %g iterations' % iter)
        break
    else:
        print('iteration: %g, L2 error = %.6f (%.6f required)' % (iter+1, error, tol))

    # Update solution at previous time step
    u0.x.array[:] = u1.x.array[:]
    k0.x.array[:] = k1.x.array[:]
    e0.x.array[:] = e1.x.array[:]

Any suggestion on how can I solve this issue is welcome :slight_smile:

P.S. the full execution of the code is also way slower then the original one (Implementing k-epsilon model in FEniCS - how to implement terms in variational formulation that depend on a solution?), any suggestion also on that?
These are the ‘petsc_options’ that I’m using for all the Linear Problems:

petsc_options = {
        "ksp_type": "cg",
        "ksp_rtol": str(tol),
        "pc_type": "hypre",
        "pc_hypre_type": "boomeramg",
        "pc_hypre_boomeramg_max_iter": 1,
        "pc_hypre_boomeramg_cycle_type": "v",
        "pc_hypre_boomeramg_print_statistics": 0, # 1 <-- to print all the informations
    }
2

Start by replacing all interative options by a direct solver, to ensure that the code is converging, i.e.

petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps", "ksp_error_if_not_converged": True}

Secondly, why do you implement all the definitions of the variational forms within the loop? I would strongly recommend defining all variational forms outside of any lopping, as you are not re-using matrices at all.

Finally, for debugging purposes I would recommend trying with a “normal” inlet condition and normal outlet, to eliminate any issues introduced by MPC.

Also, you are saying you get almost the same results, but do not actually show this by visualization, which makes it even harder to give you guidance.

1 Like
3

thanks, by using the default options now the code works properly and gives me the expected results :wink:

However, you suggested to define all variational forms outside the loop; I’ve already tried, but when I’m doing that the code seems to not work; indeed, it tells me that the error = 0 just after 1 iteration (because u1 remain equal to zero, like u0 is at the very first iteration)

I’ll put here an example of when I put the variational forms outside the loop:


# Defining the variational forms:

# N-S in 3 steps (Chorin method):
# Step 1: Tentative velocity step (u*):
F1  = (1. / dt) * dot(u - u0, v)*dx \
    + dot(dot(u0, nabla_grad(u)), v)*dx \
    + (nu + nut) * inner(grad(u), grad(v))*dx 

a_1 = form(lhs(F1)); l_1 = form(rhs(F1))


# Step 2: Pressure corrrection step (p):
a_2 = dot(grad(p), grad(q))*dx; l_2 = (-1. / dt) * dot(div(u1), q)*dx 


# Step 3: Velocity correction step (u):
a_3 = form(dot(u, v)*dx); l_3 = form(dot(u1, v)*dx - dt * dot(grad(p1), v)*dx)


# k-e turbolence model:
# Updating the variables of the model
nut, prod_k, react_k, prod_e, react_e = k_epsilon_variables()

# kinetic (k) equation:
FK  = (1. / dt) * dot(k - k0, phi)*dx \
    + dot(dot(u1, nabla_grad(k)), phi)*dx \
    + (nu + nut / 1.0) * inner(grad(k), grad(phi))*dx \
    - dot(prod_k, phi)*dx + dot(react_k * k, phi)*dx 

a_k = form(lhs(FK)); l_k = form(rhs(FK))


# epsilon () equation:
FE  = (1. / dt) * dot(e - e0, psi)*dx \
    + dot(dot(u1, nabla_grad(e)), psi)*dx \
    + (nu + nut / 1.3) * inner(grad(e), grad(psi))*dx \
    - dot(prod_e, psi)*dx + dot(react_e * e, psi)*dx 

a_e = form(lhs(FE)); l_e = form(rhs(FE))


start_case = time.time()
# main loop
iter_max = 1000 # 1000 (Remember that: dt = 0.05 (1000 * 0.05 sec = 50 sec)
tol = 1e-5

petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps", "ksp_error_if_not_converged": True}
for iter in range(iter_max):
    # N-S in 3 steps (Chorin method):
    # Step 1: Tentative velocity step (u*):
    problem1 = dolfinx_mpc.LinearProblem(a_1, l_1, mpc_u, bcs=bcu)
    u1 = problem1.solve()

    # Step 2: Pressure corrrection step (p):
    problem2 = dolfinx_mpc.LinearProblem(a_2, l_2, mpc_p, bcs=bcp)
    p1 = problem2.solve()

    # Step 3: Velocity correction step (u):
    problem3 = dolfinx_mpc.LinearProblem(a_3, l_3, mpc_u, bcs=bcu)
    u1 = problem3.solve()

    # k-e turbolence model:
    # Updating the variables of the model
    nut, prod_k, react_k, prod_e, react_e = k_epsilon_variables()

    # kinetic (k) equation:
    problem_k = dolfinx_mpc.LinearProblem(a_k, l_k, mpc_k, bcs=bck)
    k1 = problem_k.solve()

    # epsilon () equation:
    problem_e = dolfinx_mpc.LinearProblem(a_e, l_e, mpc_e, bcs=bce)
    e1 = problem_e.solve()
    # assert isinstance(e1, Function)

    # Make sure k and eps are positive
    low_bound = 1e-16
    dimension = len(k1.x.array)
    k1.x.array[:] = np.max([k1.x.array[:], low_bound * np.ones(dimension)], axis=0)  
    e1.x.array[:] = np.max([e1.x.array[:], low_bound * np.ones(dimension)], axis=0)     

    # Check if converged
    err = u0 - u1
    error_form = form(dot(err, err) * dx)
    error = assemble_scalar(error_form)
    # error = np.sqrt(assemble_scalar(error_form))
    if error < tol:
        print('Fully developed flow has been reached in %g iterations' % iter)
        break
    else:
        print('iteration: %g, L2 error = %.6f (%.6f required)' % (iter+1, error, tol))

    # Update solution at previous time step
    u0.x.array[:] = u1.x.array[:]
    u0.x.scatter_forward()
    k0.x.array[:] = k1.x.array[:]
    k0.x.scatter_forward()
    e0.x.array[:] = e1.x.array[:]
    e0.x.scatter_forward()

end_case = time.time()
print(f"\n    -> Converged in {iter} iterations (t: {end_case - start_case:.2f} sec) \n")

the output gives me:

(fenicsx-new) gianlucamongelli@Mac tubolence_modeling % python a_unit_sqaure_k_e_copy.py 
iteration: 1, L2 error = 400.000000 (0.000010 required)
Fully developed flow has been reached in 1 iterations

    -> Converged in 1 iterations (t: 0.09 sec)

Am I doing something wrong inside the loop?