Problem with PETSc solver. error code 77

Linear Algebra
1

Hello, I am running a nonlinear structural problem (Growth and remodeling of soft tissues). right now, my problem is that I am facing the error

"Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 77 (Petsc has generated inconsistent data).
*** Where:   This error was encountered inside /home/conda/feedstock_root/build_artifacts/fenics-pkgs_1617882212586/work/dolfin/dolfin/la/PETScKrylovSolver.cpp."

‘’’
I have searched a little bit, and I could not find a proper answer to this problem. There are only two topics which one of them addresses the PETSc version. I have checked my PETSc version and Fenics version, and they are up to date. Another answer was maybe it is because of lack of memory. I am using 150 gigs of memory, and when I am monitoring my memory usage, it is less than 1 percent of 150 gigs. I could upload my code for more details. But if anybody could help me to solve my problem that would be great.

Thanks

1 Like
2

Please provide a minimal code example that reproduces the error. Otherwise it is impossible to help you.

3

Thank you for your quick response. I should explain what I am doing. I am solving the dynamics of a soft tissue (continuum mechanics equation) coupled with the thickening and stiffening of the soft tissue. So I have two different cycles. one dynamic cycle which in my main loop I use a Flag_cardiac_cycle and second, Flag_local_growth which I thicken and stiffen my material.

Here is a part of my code:
weak form:


###########################################################

'''Making weak form'''


d_psi_bar_d_C_bar = ( C10Y + C0_calc*C0 * C_H * C1_H * exp(C1_H*(I1_iso - 3.0)) ) * I
d_psi_bar_d_C_bar = d_psi_bar_d_C_bar + ( C0*C0_calc * C0_H * C2_H * (1.0 - (1./(I4_iso**.5)) ) * exp(C2_H* ( (I4_iso**.5) - 1.0)**2)  ) *Dir_tensor
S_isc = 2.0*J_e**(-2.0/3.0)*Dev(d_psi_bar_d_C_bar, C_e)
second_PK_stress =  S_vol + S_isc
vel = (gamma/(beta2*dt))*(ub - ub0) - (gamma/beta2 - 1.0)*vel0 - dt*(gamma/(2.0*beta2) - 1.0)*a0 #vel and acceleration at t_(n+1)
a  = (1.0/(beta2*dt**2))*(ub - ub0 - dt*vel0) - (1.0/(2.0*beta2) - 1.0) * a0
if (Flag_generalized_alpha ): #generalized alpha integrator
 a = (1.0 - alpha_m) * a + alpha_m * a0
 vel = (1.0 - alpha_f) * vel + alpha_f * vel0
Side_flag = Constant(1.0)
Functional_b_isotropic = inner(dot(F_e,second_PK_stress),grad(vb))*dx() + Side_flag* P*J_e*inner(inv(F_e.T)*n_function,vb)*ds(1) #+P.. will apply pressure into the surface p is in line 690
Side_flag2 = Constant(0.0)
Functional_b_isotropic = Functional_b_isotropic + Side_flag2 * P*J_e* inner(inv(F_e.T)*n_function,vb)*(ds(2))#presssure on the aortic side of Leaflet
Functional_b_isotropic = Functional_b_isotropic +  P_wall*J_e*inner(inv(F_e.T)*n_function,vb)* (ds(4)) 
Functional_b_isotropic = Functional_b_isotropic +  P_wall1*J_e*inner(inv(F_e.T)*n_function,vb)* (ds(3)) 
Functional_b_isotropic = Functional_b_isotropic + rho * dot(a,vb) * dx() +damp*dot(vel,vb)*dx()
Functional_b_isotropic = Functional_b_isotropic + beta/h * dot(ub,n_function) * dot(vb,n_function)* ds(5) + beta/h * dot(ub,n_function) * dot(vb,n_function)* ds(8)
J_b_isotropic = derivative(Functional_b_isotropic, ub, du_b)
if rank == 0:
    print ('Assembling solver...')
problem_isotropic = NonlinearVariationalProblem(Functional_b_isotropic, ub, bcs_b, J_b_isotropic)
solver_isotropic = NonlinearVariationalSolver(problem_isotropic)



#################################
######my main loop: 


while cycle_number <N_cycles:
	if cycle_number==0:
		count_timestep=0
		t=0
	else:
		t=0

	vel0.vector()[:] = vel0.vector() * 0.
	a0.vector()[:] = a0.vector() * 0.
	xdmffile_u = XDMFFile(output_dir + 'displacement_'+str(cycle_number)+'.xdmf')
	xdmffile_s = XDMFFile(output_dir + 'stress_'+str(cycle_number)+'.xdmf')
	xdmffile_strain = XDMFFile(output_dir + 'strain_'+str(cycle_number)+'.xdmf')
	xdmffile_growth_rate = XDMFFile(output_dir + 'GR_'+str(cycle_number)+'.xdmf')
	xdmffile_strain_ave = XDMFFile(output_dir + 'ave_strain_'+str(cycle_number)+'.xdmf')
	xdmffile_C0_glob = XDMFFile(output_dir + 'c0glob_'+str(cycle_number)+'.xdmf')
	xdmffile_c0_loc = XDMFFile(output_dir + 'c0_loc_'+str(cycle_number)+'.xdmf')
  ###########################################################
  ###cardiac cycle loop
	if cycle_number%2==0:
		Flag_global_growth=False
		Flag_local_growth=False
		Flag_cardiac_cycle=True
		count_timestep=0
		temp_strain_calc=np.zeros((element_numbers))
	else:
		Flag_global_growth=False
		Flag_local_growth=True
		Flag_cardiac_cycle=False

  ###########################################################
	if (Flag_cardiac_cycle):
		rho.assign(1100)
		while t<=0.42:
			if(1):
				if t<0.0:
					dt_s=dt_coarser
				elif t<0.04:
					dt_s=dt_med
				elif t< 0.25:
					dt_s=dt_med
				elif t<=0.419:
					dt_s=dt_coarse
				else:
					dt_s=dt_fine
			t_prev=t
			t+=dt_s
			dt.assign(dt_s)
			temp_array_P = P.vector().get_local()
			temp_array_P_wall = P_wall.vector().get_local()
			temp_array_P_wall1 = P_wall1.vector().get_local()
			T_alpha = (1.0 - alpha_f) * t + (alpha_f * t_prev)
			P_current_trans = np.interp(T_alpha,Time_array,P_trans_array)
			P_current_wall = np.interp(T_alpha,Time_array,P_aort_array)
			P_current_wall1 = np.interp(T_alpha,Time_array,P_vent_array)
			shift=np.interp(T_alpha,Time_array2,P_shift)
			if (P_current_trans < 0 ):
				Side_flag.assign(0.0)
				Side_flag2.assign(1.0)
				P_current_trans_in = -P_current_trans
				vent_side=False
			else:
				P_current_trans_in=P_current_trans#+cycle_number*shift
				Side_flag.assign(1.0)
				Side_flag2.assign(0.0)
				vent_side=True
			if rank==0:
				print('cardiac_cycle')
				print ('cycle=', cycle_number)
				print ('time=', t)
				print ("tot_pressure", P_current_trans_in)
				if vent_side:
					print ("pressure_in_vent=", P_current_trans)
					print ("shifting_pressure", shift)
				else:
					print ("pressure_in aort=", P_current_trans)
			temp_array_P[:] = P_current_trans_in
			temp_array_P_wall1[:] = P_current_wall1
			temp_array_P_wall[:] = P_current_wall
			P.vector().set_local(temp_array_P)
			P_wall.vector().set_local(temp_array_P_wall)
			P_wall1.vector().set_local(temp_array_P_wall1)
			solver_isotropic.solve()
			update(ub,ub0,vel0,a0,beta2,gamma,dt_s)
			solve(a_m2 == L_m2, strain_eps_proj)
			temp_strain_calc+=strain_eps_proj.vector().get_local()
			# temp_strain_calc=strain_eps_proj.vector().get_local()
			if count_saving_file1%7==0:
				ub.rename("displacement","Label")
				xdmffile_u.write(ub, t)
				strain_eps_proj.rename("strain","Label")
				xdmffile_strain.write(strain_eps_proj,t)
				stress_tensor = (1./J_e)* F_e*second_PK_stress*F_e.T
				TTTT=TensorFunctionSpace(bmesh, 'DG', 0)
				eqSigma = Function(TTTT)
				eqSigma.assign(project(stress_tensor,TTTT))
				eqSigma.rename("stress","Label")
				xdmffile_s.write(eqSigma, t)
			count_saving_file1+=1
			count_timestep+=1
		temp_strain_calc=temp_strain_calc/count_timestep
		for kkk in range(count_f):
			if abs(temp_strain_calc[Fibrosa_IDs[kkk]])>strain_threshhold:
				calcif_tags[Fibrosa_IDs[kkk]]=1
				# cell_number_calc = Cell(bmesh, Fibrosa_IDs[kkk])
				point_1=(V0.dofmap().cell_dofs(Fibrosa_IDs[kkk])[0])
				point_2=(V0.dofmap().cell_dofs(Fibrosa_IDs[kkk])[1])
				point_3=(V0.dofmap().cell_dofs(Fibrosa_IDs[kkk])[2])
				point_4=(V0.dofmap().cell_dofs(Fibrosa_IDs[kkk])[3])
				growth_local_rate_numpy[point_1]=1
				growth_local_rate_numpy[point_2]=1
				growth_local_rate_numpy[point_3]=1
				growth_local_rate_numpy[point_4]=1
		growth_rate_local.vector().set_local(growth_local_rate_numpy)
		growth_rate_local.vector().apply('insert')
		growth_rate_local.rename("gr","Label")
		xdmffile_growth_rate.write(growth_rate_local)
	elif (Flag_local_growth):
		loc_flag.assign(1.0)
		rho.assign(0.0)
		C_0_temp_calc[:]=C0_calc.vector().get_local()[:]
		for hh in range (element_numbers):
			if calcif_tags[hh]==1:
				if C_0_temp_calc[hh]<=2:
					C_0_temp_calc[hh]=C_0_temp_calc[hh]*stiffness_factore_local
		C0_calc.vector().set_local(C_0_temp_calc[:])
		C0_calc.rename("c0loc","Label")
		xdmffile_c0_loc.write(C0_calc)
		for GG in range (601):
			if rank==0:
				print("loc_growth_cycle", GG)
			growth_constant_loc.assign(GG*.001)
			solver_isotropic.solve()
			if GG%50==0:
				ub.rename("displacement","Label")
				xdmffile_u.write(ub, GG)
		
		
	cycle_number=cycle_number+1
	if rank ==0:
		print("cycle_number",cycle_number)

#########################

here the number of cycles is 24. the code works without any problem in cycle_number 0,1 and 2 (no matter if we are in Flag_local_growth or Flag_cardiac_cycle loop. However, when it goes to cycle 3 ( Flag_local_growth become true), I am getting the error I mentioned. I have tried to print the linear solver residual and this is the last line of the residual;

110 KSP preconditioned resid norm 4.341064278800e-09 true resid norm 1.159734350726e-05 ||r(i)||/||b|| 2.784807915884e-05

This is my preconditioner for my solver:

prm = solver_isotropic.parameters
#print prm["newton_solver"]["maximum_iterations"]
prm["newton_solver"]["maximum_iterations"] = 200
prm["newton_solver"]["krylov_solver"]["nonzero_initial_guess"] = True
prm['newton_solver']['absolute_tolerance'] = 1E-6 
prm['newton_solver']['linear_solver'] = 'tfqmr' #'gmres' bicgstab
prm['newton_solver']['preconditioner'] = 'petsc_amg'
prm["newton_solver"]["krylov_solver"]["monitor_convergence"] = True
prm['newton_solver']['krylov_solver']['maximum_iterations'] = 200000

If you have more other questions, just let me know.

4

When I asked for a minimal code example, it has to be a code that I can execute on my system, that produces the error you are experiencing.

Just having a piece of the variational form is not sufficient.

There are a list of examples of good minimal code examples:

Please also write which version of FEniCS and PETSc you are using.

5

I can upload the whole code but because my geometry and boundary conditions and conditions I have for my problem are complex, I could not make it a minimal code example. Also, my main code takes about 900 minutes in 64 processors to reach the point that I get errors. The conditions I have it is somehow impossible to implement it in a simple model. Actually, I have solved my problem with simpler conditions and geometry. If you don’t mind and have extra time I can send you my geometry and the main code. I need your email if it is possible for you.

6

Well, so the boundary conditions should not be dependent if the number of the cells, so you should be able to use a coarser geometry, with fewer cells. which in turn will be faster to simulate to check that the code actually runs.

It is impossible for anyone to debug or help you if the code requires 900 minutes of simulation to get to an error.

You need to try to simply your code, step by step, and see if you can find the cause of your error.

For instance, at what time-step does the problem fail. Do you complete a certain number of time steps? How does your solution look at the time step before failure? Are any new terms or integration domains introduced between the second to last and last time step?

1 Like
7

I know what you mean. because of that, I have just asked for the code number error (PETSc error code is: 77 (Petsc has generated inconsistent data). my simulation has two loops and individually, each loop works fine when I am running them individually (simpler condition). In my final version of the code, I need to couple these two loops, (calculating strain from the first loop, and update my geometry and material properties based on that strain in the second loop, and again calculating strain from the updated material properties, and so on). So, every other iteration it goes through each loop. I don’t have a problem with the 1st 2nd and 3rd iteration, and my problem starts at the beginning of the 4th iteration when my code enters in the second loop for the second time (updating material properties from the strain). I don’t have time integration in that loop. and by printing the residuals for Newtonian and linear solvers I can see my code is converging. however, suddenly it gives me the error. I have checked my results till the last point before this error comes up and everything looks fine. Anyway if you have any other suggestions, I will be happy to hear from you.

8

Out of curiosity, how many cells are there in your mesh, and how many degrees of freedom in your problem?

9

I have 323930 cells and 64922 points.

10

I wouldn’t call that a massive mesh. It should be possible to run a problem with that amount of cells on a laptop/desktop with 1-5 MPI processes.

It would help if you could clarify how many degrees of freedom, what spaces you are using, and how many variational forms (including projection and interpolation) that you are using.

11

ok, let me tell you, I am calculating strain element wise and it is just a one direction of strain so it will be a FunctionSpace(mymesh,‘DG’,0), and for each element if my strain value is more than a threshhold, I am tagging that element. then I should change one of the material properties which is based on same space (FunctionSpace(mymesh,‘DG’,0)) for that element, while another property I should calculate them node wise for that tagged elements (4 points in each element) which is based on the (VectorFunctionspace (mymesh,'CG,1)). So if you look into my part of code that it is end of the cardiac cycle loop, I have tagged two different parameters (one node wise, one element wise).

for kkk in range(count_f):
			if abs(temp_strain_calc[Fibrosa_IDs[kkk]])>strain_threshhold:
				calcif_tags[Fibrosa_IDs[kkk]]=1  ####### this is a material property that I have changed element wise
				# cell_number_calc = Cell(bmesh, Fibrosa_IDs[kkk]) 
				point_1=(V0.dofmap().cell_dofs(Fibrosa_IDs[kkk])[0]) ##### here I have tagged each element points for second material property which is based on the vectorfunctionspace
				point_2=(V0.dofmap().cell_dofs(Fibrosa_IDs[kkk])[1])
				point_3=(V0.dofmap().cell_dofs(Fibrosa_IDs[kkk])[2])
				point_4=(V0.dofmap().cell_dofs(Fibrosa_IDs[kkk])[3])
				growth_local_rate_numpy[point_1]=1
				growth_local_rate_numpy[point_2]=1
				growth_local_rate_numpy[point_3]=1
				growth_local_rate_numpy[point_4]=1
		growth_rate_local.vector().set_local(growth_local_rate_numpy)
		growth_rate_local.vector().apply('insert')
		growth_rate_local.rename("gr","Label")
		xdmffile_growth_rate.write(growth_rate_local)

Then I update them in my second loop for the both of them:

for hh in range (element_numbers):
			if calcif_tags[hh]==1:
				if C_0_temp_calc[hh]<=2:
					C_0_temp_calc[hh]=C_0_temp_calc[hh]*stiffness_factore_local
		C0_calc.vector().set_local(C_0_temp_calc[:])
		C0_calc.rename("c0loc","Label")

for GG in range (601):
			if rank==0:
				print("loc_growth_cycle", GG)
			growth_constant_loc.assign(GG*.001)
			solver_isotropic.solve()

So then I could not get the problem solved in the last part solver_isotropic.solve()

I hope I could explain it clearly.

12

If you are doing things element wise, i.e. looping over facets and cells, I would strongly suggest you using the C++ interface of dolfin as your mesh has 300 000 cells. Doing such loops and individual solves in Python is alot slower than doing them in C++

1 Like
13

Thank you for your great information, I found out that my code diverges in Newtonian iteration solver. do you have any suggestions to prevent that? generally, what do people do to help the convergence of Newtonian solver?

14

@nate has written some good tips at: Default absolute tolerance and relative tolerance - #4 by nate

2 Likes