distutils.errors.DistutilsExecError

dolfinx
1

Hello all,

I am working on a phase field problem. The formulation I use is :
image
[ref : Maurini and Farell ]

I got pretty good results with w(alpha) = alpha² but with w(alpha) = alpha, I met this issue :


libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c: In function 'functionspace_form_libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746_0':
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1173:1: warning: parameter names (without types) in function declaration
 1173 | static ufcx_function_space functionspace_Max(psi_plus) =
      | ^~~~~~
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1173:28: error: invalid storage class for function 'functionspace_Max'
 1173 | static ufcx_function_space functionspace_Max(psi_plus) =
      |                            ^~~~~~~~~~~~~~~~~
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1175:1: error: field name not in record or union initializer
 1175 | .finite_element = &element_703c39b3216cbd499bb1267d13b485382ccc3f67,
      | ^
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1175:1: note: (near initialization for '(anonymous)')
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1176:1: error: field name not in record or union initializer
 1176 | .dofmap = &dofmap_703c39b3216cbd499bb1267d13b485382ccc3f67,
      | ^
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1176:1: note: (near initialization for '(anonymous)')
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1177:1: error: field name not in record or union initializer
 1177 | .geometry_family = "Q",
      | ^
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1177:1: note: (near initialization for '(anonymous)')
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1178:1: error: field name not in record or union initializer
 1178 | .geometry_degree = 1,
      | ^
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1178:1: note: (near initialization for '(anonymous)')
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1179:1: error: field name not in record or union initializer
 1179 | .geometry_basix_cell = 4,
      | ^
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1179:1: note: (near initialization for '(anonymous)')
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1180:1: error: field name not in record or union initializer
 1180 | .geometry_basix_variant = 0
      | ^
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1180:1: note: (near initialization for '(anonymous)')
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1187:11: warning: implicit declaration of function 'functionspace_Max' [-Wimplicit-function-declaration]
 1187 |   return &functionspace_Max(psi_plus);
      |           ^~~~~~~~~~~~~~~~~
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1187:29: error: 'psi_plus' undeclared (first use in this function)
 1187 |   return &functionspace_Max(psi_plus);
      |                             ^~~~~~~~
libffcx_forms_d794fd15c146f0ff61fe6d1ac46d700b634b8746.c:1187:29: note: each undeclared identifier is reported only once for each function it appears in
Traceback (most recent call last):
  File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/setuptools/_distutils/unixccompiler.py", line 117, in _compile
    self.spawn(compiler_so + cc_args + [src, '-o', obj] +
  File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/setuptools/_distutils/ccompiler.py", line 917, in spawn
    spawn(cmd, dry_run=self.dry_run, **kwargs)
  File "/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/lib/python3.10/site-packages/setuptools/_distutils/spawn.py", line 68, in spawn
    raise DistutilsExecError(
distutils.errors.DistutilsExecError: command '/stck/lmersel/CODE_PROGRAMME/anaconda3/envs/fenics/bin/x86_64-conda-linux-gnu-cc' failed with exit code 1

Does that related to this discussion here ? But the difference is that I am not using mixed space like the discussion mentioned and also the code worked on Fenics 2019.

My variational formulation on Fenicsx is

#---------------- Mechanical equation
residual = ufl.inner( sigma_amor(u, alpha, G , K0, ndim) , epsilon(v) ) * dx# - dot( tn, v ) * ds #residual
jacobian = ufl.derivative(residual, u, u_trial) # jacobian

#--------------- Phase Field equation
f_alpha = H_new*ufl.inner(beta, gk_prime(alpha)) * dx \
          + (Gc/(cw*lw)) * (beta*w_prime(alpha) + (2.*lw**2)*ufl.inner(ufl.grad(beta), ufl.grad(alpha))) * dx 
J_alpha = ufl.derivative(f_alpha, alpha, alpha_trial) # jacobian

#---------------- RESOLUTION
t= 0.
tol = 0.001 # 1e-6
maxiter = 300
begin_tt = 1
it_time = begin_tt

for t in time_step[begin_tt:2]:
    
    u_app.t = t
    u_app_func.interpolate(u_app)
    loading[it_time]= np.array([it_time,u_app.t])
    print("\n================================================================")
    print("\n -> Time increment: %d, Time: %.3g s, Ux_given: %.3g mm  "%(it_time, t, u_app.t))
    print("\n================================================================")

    it_staggered = 0
    res_H = 1.
    err = 1.
    print("\n ----- > Start Staggered Loop")
    while (res_H > tol) and (it_staggered < maxiter):
        it_staggered +=1

        #----Step 1 : Displacement resolution
        print("\n NRSolver for the displacement pb ")
        print("_____________________________________")

        bc_left = fem.dirichletbc( fem.Constant(msh, ScalarType(0.)), facet_left_dofs,V_u.sub(1))
        bc_right = fem.dirichletbc( fem.Constant(msh, ScalarType(0.)), facet_right_dofs,V_u.sub(1))
        bc_bottom = fem.dirichletbc( u_zero,facet_bottom_dofs, V_u)
        bc_top = fem.dirichletbc( u_app_func, facet_top_dofs)
        bcu = [bc_top, bc_right, bc_bottom,bc_left ]

        #Build homogeneous bc
        u_app_0 = TopLoad(t=0.)
        u_app_func_h.interpolate(u_app_0)
        bc_top_h = fem.dirichletbc( u_app_func_h, facet_top_dofs)
        bcu_h = [bc_top_h, bc_right, bc_bottom,bc_left]

        #Newton Raphson resolution
        iter_u = newton_solver( u,  res = residual, bcs = bcu,bcs_hom = bcu_h, jac = jacobian)
        
        #----Step 2 : Computation of H energy history

        Pel_pos_expr = fem.Expression(psi_plus(u,K0,G), V_H.element.interpolation_points)
        psi_elastic_pos.interpolate(Pel_pos_expr)
        H_new.x.array[:] = np.max((psi_elastic_pos.x.array[:], H_inc.x.array[:]),axis=0)

        

        #----Step 3 : Damge resolution

        print("\n NRSolver for damage pb")
        print("_____________________________________")
        
        bcalpha = [fem.dirichletbc( fem.Constant(msh, PETSc.ScalarType(1.)), facet_notch_dofs,V_alpha)]

        #Build homogeneous bc
        bcalpha_h = [fem.dirichletbc( fem.Constant(msh, PETSc.ScalarType(0.)), facet_notch_dofs,V_alpha)]

        #Newton Raphson resolution
        iter_alpha = newton_solver( alpha,  res = f_alpha, bcs = bcalpha,bcs_hom = bcalpha_h, jac =  J_alpha)

        print("_____________________________________")

        H_gap = H_new.vector.array - H_old.vector.array
        index_max,H_gap_max = np.argmax(H_gap), H_gap.max()
        res_H = abs(H_gap_max)/ H_old.vector.array[index_max]
        
        print("Staggered Iter:", it_staggered)
        print("abs(H_gap_max)/ H_old.vector.array[index_max]", res_H)

        #---- Step 4 : Update previous solution
        u_old.x.array[:] = u.x.array[:]
        alpha_old.x.array[:] = alpha.x.array[:]
        g = gk(alpha_old,1.0e-6)
        alpha_old.x.array[:] = alpha.x.array[:]
        H_old.x.array[:] = H_new.x.array[:]
 
    print("\n ----- > END Staggered Loop")
    H_inc.x.array[:] = H_new.x.array[:]

    # monitor the results    
    if comm.rank == 0:
        print( "\nalpha_max:", (alpha.x.array.max() ))
        print("\nIterationU and IterationAlpha: (%i,%i)"%(iter_u,iter_alpha))
        print("\nElastic and surface energies : (%g,%g)"%(elastic_energy_scalar,surface_energy_scalar - surf_off_set)) 


    it_time +=1

And the newton-raphson solver is :

def newton_solver(u, res, bcs,bcs_hom , jac, max_iter=100, rtol=1e-10, atol=1e-6 ):

    #applied BC + homogenize BCs
    fem.petsc.set_bc(u.vector, bcs)
    bcs = bcs_hom

    jac_form, res_form = fem.form(jac), fem.form(-res)
    J = fem.petsc.assemble_matrix(jac_form, bcs = bcs)
    J.assemble()
    b = fem.petsc.create_vector(res_form)
    delta_u = J.createVecRight()

    for k in range(max_iter):

        # Assemble the system for the Newton update with boundary conditions

        with b.localForm() as loc_1:
            loc_1.set(0)
        b = fem.petsc.assemble_vector(res_form)
        fem.petsc.apply_lifting(b, [jac_form], [bcs]) #,x0=u.vector, scale=-1.0) 
        b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
        fem.petsc.set_bc(b, bcs)#,x0=u.vector, scale=-1.0) 

        # Create CG Krylov solver and turn convergence monitoring on
        solver = PETSc.KSP().create(MPI.COMM_WORLD)
        pc = solver.getPC()
        solver.setType('preonly')
        pc.setType('lu')
        solver.setConvergenceHistory()

        # Set matrix operator
        solver.setOperators(J)

        # Compute solution
        solver.setMonitor(lambda ksp, its, rnorm: print("Iteration: {}, rel. residual: {}".format(its, rnorm)))
        solver.solve(b, delta_u)

        # Check residual n°1
        residual = J * delta_u - b
        print(f"The relative residual is: {residual.norm() / b.norm()}.")

        # Update the solution
        u.x.array[:] = u.x.array[:] + delta_u.array[:]

        print("norm(b) : ", b.norm())
        
        if   residual.norm() /b.norm() < rtol and k < max_iter:
            break


    else:
        raise RuntimeError("could not converge, norm_u {}, norm_delta_u {}".format(norm_u, norm_delta_u))
    
    
    return k

Thank you for the help !

Best regards,

Lamia

2

Please rewrite the code in your post such that it can reproduce the error message.

This means:

  • Creating the smallest possible lines of code that can be copy-pasted and executed by anyone else, and reproduce your error message.

To do so, I strongly suggest using any of the built-in meshes (create_unit_square/create_unit_cube) etc. and remove all lines of code after the error message appears.

This means that you should not need any solver declarations, as the problem happens in the form compiler step (or so it seems).

3

It’s really surprising, the error has completetly disappeared. And I have done nothing less than running the same code again. It seems that it fixed but it is a bit frustrating because I don’t know how.
But thank you for the answer ! I will try to follow the guidelines for the future posts.