Diverging results with finer meshes @ magneto-mechanical, periodical, nonlinear simulation

1

Hello,

I do some coupled magneto-mechanical, nonlinear, 3D simulation on a cubical periodic domain of size L^3. I am currently facing this problem: For coarse meshes (up to 55k elements) the solution is fine and meets the expectations. However, contain the meshes more elements (>75k) the solution somehow gets on a wrong path and finally diverges (see figure).

PMech-H_Fenics
PMech-H_Fenics.png1120×840 48.4 KB
This problem does not appear if a purely mechanical or magnetic analysis is run, neither if the domain is not periodically implemented. Then the solutions are similar, no matter how many elements.

I use tetrahedrons as elements. A scalar potential approach is used to compute the magnetic quantities.

The periodic boundary conditions are this way implemented and seem to be working fine:

class PeriodicDomain(SubDomain):	
	# return true for left, bottom and front (x = 0, y = 0, z = 0)
	# but not edges shared with sides to map
	def inside(self, x, on_boundary):
		# on left, bottom or front side
		bool_1 = bool( (x[0] < tol) or (x[1] < tol) or (x[2] < tol) )
		# on shared edges of x=0-side
		bool_2 = bool( ((x[0] < tol) and (x[2] > L - tol)) or ((x[0] < tol) and (x[1] > L - tol)) )
		# on shared edges of y=0-side
		bool_3 = bool( ((x[1] < tol) and (x[0] > L - tol)) or ((x[1] < tol) and (x[2] > L - tol)) )
		# on shared edges of z=0-side
		bool_4 = bool( ((x[2] < tol) and (x[0] > L - tol)) or ((x[2] < tol) and (x[1] > L - tol)) )
		# combine bools
		return bool_1 and not(bool_2) and not(bool_3) and not(bool_4) and on_boundary
		
	# map right, top, back side to sides defined by inside
	# map edges shared by this sides extra
	def map(self, x, y):
		# edge between z = L and x = L
		if (x[0] > L - tol) and (x[1] > L - tol) and (x[2] > L - tol):
			y[0] = 0
			y[1] = 0
			y[2] = 0
		elif (x[0] > L - tol) and (x[2] > L - tol):
			y[0] = 0
			y[1] = x[1]
			y[2] = 0
			# y[0] = x[0] - L
			# y[1] = x[1]
			# y[2] = x[2] - L
			cnt[0] += 1
		# edge between x = L and y = L
		elif (x[0] > L - tol) and (x[1] > L - tol):
			# y[0] = x[0] - L
			# y[1] = x[1] - L
			# y[2] = x[2]
			y[0] = 0
			y[1] = 0
			y[2] = x[2]
			cnt[1] += 1
		# edge between y = L and z = L
		elif (x[1] > L - tol) and (x[2] > L - tol):
			# y[0] = x[0]
			# y[1] = x[1] - L
			# y[2] = x[2] - L
			y[0] = x[0]
			y[1] = 0
			y[2] = 0
			cnt[2] += 1
		# right side x = L
		elif x[0] > L - tol:
			# y[0] = x[0] - L
			# y[1] = x[1]
			# y[2] = x[2]
			y[0] = 0
			y[1] = x[1]
			y[2] = x[2]
			cnt[3] += 1
		# top side y = L
		elif x[1] > L - tol:
			# y[0] = x[0]
			# y[1] = x[1] - L
			# y[2] = x[2]
			y[0] = x[0]
			y[1] = 0
			y[2] = x[2]
			cnt[4] += 1
		# back side z = L
		elif x[2] > L - tol:
			# y[0] = x[0] 
			# y[1] = x[1]
			# y[2] = x[2] - L
			y[0] = x[0] 
			y[1] = x[1]
			y[2] = 0
			cnt[5] += 1
		# should not happen, catch case and prompt to user
		else:
			y[0] = 2000
			y[1] = 0
			y[2] = 0
			cnt[6] += 1
			print("\n!!!\n>>>>>Error within assignment of periodic boundaries.\n>>>>>Consider using different tolerance.\n!!!\n")
pbc	= PeriodicDomain()

The function spaces and elements are implemented this way:

V_fluc = VectorElement('P', tetrahedron, 2)                 # vector elements for fluctuation field
V_pot  = FiniteElement('P', tetrahedron, 2)                 # (scalar) finite elements for magnetic potential
FE     = MixedElement([V_fluc, V_pot])                      # combine elements to mixed finite elements
W      = FunctionSpace(mesh, FE, constrained_domain = pbc)

Additionally one node is fixed in displacement and the potential set to 0. I use the direct ‘mumps’ solver with a relative tolerance of 1E-6. Fenics is used in the version 2018.1.0