Vectorfunction space degree and Quadrature degree


228
views
0
5 months ago by
Ovais  
Hi there

Is it always necessary for the quadrature degree to be equal to be VectorFunctionSpace degree ?

I am trying to solve a nonlinear elasticity problem which gets solved when i set VectorFunctionSpace degree = 2 and Quadrature degree = 4 but the solution process gets killed when i used the degree of both VectorFunctionSpace as well as Quadrature = 4.

Have not posted code because that makes the question very long.

Thanks. (please also comment if that is not required or necessary as a rule)
Community: FEniCS Project
1
Can you post a minimum working example? Not sure what you mean with quadrature degree here.
written 5 months ago by Miguel  

mesh = Mesh('biv_mesh.xml')

def Make_SubDomains(mesh):
	mf = FacetFunction("size_t", mesh)						
	mf.set_all(0)												
	File("facet_function.xml")>> mf
	return mf

# ----------------------------------------------------------------------------- #

import numpy as np 
import matplotlib.pyplot as plt

# Active tension assumed gaussian function with values in Pa total 901 entries
active_tension = np.loadtxt('active_tension_900ms.txt')*0.76                           
# Array from 0 to 0.9 with time spacing equal to 0.9/len(active_tension)
t = np.linspace(0,0.9,len(active_tension)) 					       
# Empty vector of 901 entries

pressure = np.zeros(len(t))								       
pressure_rv = np.zeros(len(t))

pressure[:500] = (t[:500]/0.5)*1500 				          
pressure_rv[:500] = (t[:500]/0.5)*1000

# Step step size of 0.1 
dt = Constant(t[1]) 								        

# Set up mesh boundaries and normal 
boundaries = Make_SubDomains(mesh)                # Identify boundaries as facets 
ds = ds[boundaries] 										  
n_mesh = FacetNormal(mesh)

# Space and functions 
U = VectorFunctionSpace(mesh,"Lagrange", 2) 	  # Function Space for displacement 
V = VectorFunctionSpace(mesh,"Lagrange", 1)       # Func Space for velocity				
Q = FunctionSpace(mesh,"Lagrange",1)  	          # Function Space for pressure
W = MixedFunctionSpace([U,V,Q])								   

uvp = Function(W)									   
u, v, p = split(uvp)
du, dv, dp = TestFunctions(W)

uvp = Function(W)
uvp_1 = Function(W)
u, v, p = split(uvp)
u_1, v_1, p_1 = split(uvp_1)
du, dv, dp = TestFunctions(W)

gamma_1 = Constant(22.6)
gamma_2 = Constant(0.01)

# Kinematics
d = u.geometric_dimension()	# geometric dimension of the vector field 
I = Identity(d)             				
F = I + grad(u)			# Deformation gradient     
C = F.T*F                   	# Right Cauchy-Green tensor
B = F*F.T

# Time-rates of deformation measures
F_dot = grad((u - u_1)/dt)
l = F_dot*inv(F)
d = 0.5*(l+l.T)
C_dot = 2*F.T*d*F
B_dot = l*B + B*l.T

# Invariants of deformation tensors
I_1 = tr(C)
I_1_dot = tr(C_dot)
J = det(F)

# Elasticity/material parameters
eta = 0.1  						
rho = Constant(1000.0) 			
alpha = 0.1 					
a = 2280.0*alpha 			
a_f = 1168.5*alpha 		 	
b = 9.726*0.8		 	
b_f = 15.779*0.75 		 	


FiberSpace = VectorFunctionSpace(mesh,'Quadrature',4)  
													  
f_0 = Function(FiberSpace) 	# fiber direction
s_0 = Function(FiberSpace) 	# sheet vector
n_0 = Function(FiberSpace) 	# sheet normal
Load('fiber_fd.h5','f_0',f_0)
Load('fiber_fd.h5','s_0',s_0)
Load('fiber_fd.h5','n_0',n_0)

f = F*f_0 
s = F*s_0 
n = F*n_0 

# Invariants
I_4f = inner(C*f_0,f_0)
I_4f_dot = inner(C_dot*f_0,f_0)

# Initial conditions and boundary conditions
zero_displacement = Expression(("0.0", "0.0", "0.0"))
bcr = DirichletBC(W.sub(0), zero_displacement, boundaries, 10)			
bcs = [bcr]
p0 = Constant(0); T_a = Constant(0) ; p1 = Constant(0)									

# Stress relations


sigma_elastic = a*exp(b*(I_1 - 3))*B + 2*a_f*(I_4f-1)*exp(b_f*pow(I_4f - 1, 2))*outer(f,f) - p*I   

sigma_viscous = gamma_1*exp(gamma_2*I_1_dot)*B_dot 												    

sigma_active = T_a*(outer(f,f)+eta*outer(s,s)+eta*outer(n,n))									    

# Variational formulation (Weak form)
P = J*(sigma_active+sigma_elastic+sigma_viscous)*inv(F.T)										    
eq1 = inner(P,grad(du))*dx + (rho/dt)*(inner(v - v_1 ,du)*dx) 									    
eq2 = (1/dt)*(inner(u - u_1,dv)*dx) - inner(v,dv)*dx
eq3 = inner(J-1,dp)*dx																			    
eq = eq1 + eq2 + eq3 + dot(J*inv(F).T*n_mesh*p0, du)*ds(30) + dot(J*inv(F).T*n_mesh*p1, du)*ds(20)			


Jac = derivative(eq, uvp)

# solve the problem for displacement , velocity and pressure
i=0
T = 400
while i<=T:
     print 'Now solving time-step:', i, '/', 400
     solve(eq==0,uvp,bcs,J=Jac,solver_parameters={'newton_solver':
     {'linear_solver':'mumps'}})
     u,v,p = uvp.split(deepcopy=True)
     file << u
     uvp_1.assign(uvp)
     i = i+1
	


written 5 months ago by Ovais  
I have reduced some of the code but i had to add this much so that you can understand. I just wanted to ask about the interdependence of VectorFunctionSpace and Quadrature degree. The code seems to get killed when degree of VectorFunctionSpace is increased from 2 to 4. Any idea of what happens behind the scene when we increase the degree ? ... and would be really helpful if you can suggest the ideal selection for this case. Thanks for your interest
written 5 months ago by Ovais  
My guess is that by increasing the degree of the FunctionSpace you generate higher order terms, which will then not be approximated sufficiently well with quadratures of degree 4. What happens if you increase the quadrature degree?
written 5 months ago by Lukas O.  
When i increase the degree of VectorFunctionSpace from 2 to 4 and keep Quadrature degree as 4, the code doesnt runs and gets killed. Where as for VectorFuncationSpace degree = 2 and Quadrature degree =4 the code runs and solves the problem. What does killed indicate ? does it indicate an error or does it indicate excessive computations. I am using Core i 7 with 16 GB RAM.

I have also added the code so that you and others interested in this post may have an idea of the complete problem.
Thanks for response.
written 5 months ago by Ovais  
The "killed" message is out-of-memory behaviour. You could try higher degree quadrature on a smaller mesh and compare the errors to choose a fitting quadrature degree. Switching from the direct mumps solver to something iterative like gmres could also help.
written 5 months ago by klunkean  

Thanks for the response.  I was able to find the answer to my querry by trying vector function space degree =1, 2 against Quadrature degree = 2 one by one and then compared the results.
For vector function space degree lesser than quadrature degree i was getting wrong or erroneous displacements where as  for vector function space degree = Quadrature degree the displacement that i got was correct. I verified the solution by comparing the fem solution with solution of scipy optimize (on simplified problem). Hence my finding is that Vector function space degree for Nonlinear hyperelasticity problems related to stresses should be equal to quadrature degree. I am open to criticism on my results. If someone from community has an objection or better explanation then it would definitely help me improve my conceptual understanding too.

Now..
I have been trying figure out way to somehow use mumps for vectorfunction space degree = 4 because mumps seems to be performing better than iterative solver for vectorfunction space degree =2. so far i am using the following :-

solve(eq==0 , w , bcs , J=Jac, solver_parameters={'newton_solver':{'linear_solver':'mumps'}})

​
 Can anyone suggest how to improve the performance of mumps ?. There are 2,44,000 dof in the system of equation that i am trying to solve.
Best regards.
written 5 months ago by Ovais  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »