### Matrix assembly takes too long...

32
views
0
10 days ago by

Hello,
I am solving the 2d Stokes equations with a viscoelastic constitutive equation in a square domain with some Dirichlet conditions (I use P2/P1/P1/P1 elements for the velocities/pressure/velocity gradient projections/stresses). In order to achieve numerical stability I use the log conformation representation of the stress tensor. In this formulation, one needs to solve for the logarithm of the stress tensor instead solving for the stress tensor itself. To do so, we diagonalize the stress tensor, and we also calculate two other tensors based on dot products of the velocity gradient tensor and a tensor that contains the eigenvectors of the stress tensor. Note that all these processes are implicit and all these tensors are used in the nonlinear variational form. In addition, these operations have to be performed at the Gauss points (not at the nodes of the triangles). I have written 3 functions to perform the diagonalization process, but when I use them the matrix assembly process becomes extremely slow... More specifically the matrix assembly process becomes approximately 30 times slower when I use these functions (log representation) compared to the time that the matrix assembly process takes when I solve directly for the stress tensor (and thus, I do not use these functions..).
Any idea how to speed up my code??? (I use optimization flags and mpi parallel)
Here are the diagonalization functions:
from dolfin import *

#---------------------------------------------------------------------------
# Function eigenspectum
#---------------------------------------------------------------------------

def eigenspectum( Sxx, Syx, Syy ):

# Define trace of the tensor
TrS = Sxx + Syy

# Define discriminant of the characteristic equation
Dis = (Sxx-Syy)*(Sxx-Syy) + 4.0*Syx*Syx

# Calculate eigenvalues
l1 = (TrS+sqrt(Dis))/2.0
l2 = l1-sqrt(Dis)

# Calculate eigenvectors
v1x = conditional(gt(abs(Syx), 1.e-9), Sxx-Syy+sqrt(Dis), 1.0)
v1y = conditional(gt(abs(Syx), 1.e-9), 2.0*Syx, 0.0)

# Normalize eigenvectors
V1_m = sqrt(v1x*v1x+v1y*v1y)

v1x = v1x/V1_m
v1y = v1y/V1_m

return (l1, l2, v1x, v1y)

#---------------------------------------------------------------------------
# Function stretch_rotational_tensors
#---------------------------------------------------------------------------

def stretch_rotational_tensors( Sxx, Syx, Syy, Gxx, Gyx, Gxy, Gyy ):

# Calculate eigenvalues & eigenvectors of S tensor
(l1, l2, v1x, v1y) = eigenspectum( Sxx, Syx, Syy )

# Define R & R transpose tensors
R = as_tensor([\
[v1x, -v1y],\
[v1y, v1x]\
])

RT = transpose(R)

# Define velocity gradient transpose tensor
GUT = as_tensor([\
[Gxx, Gyx],\
[Gxy, Gyy]\
])

# Define M tensor and its diagonal
M = dot(dot(RT,GUT),R)

Md = as_tensor([\
[M[0,0], 0.0 ],\
[0.0, M[1,1]]\
])

# Define B tensor
B = dot(dot(R,Md),RT)

# Define B tensor
Bxx = conditional(gt(abs(l1-l2), 1.e-9), B[0,0], Gxx )
Byx = conditional(gt(abs(l1-l2), 1.e-9), B[0,1], 0.5*(Gyx+Gxy))
Byy = conditional(gt(abs(l1-l2), 1.e-9), B[1,1], Gyy )

# Define W tensor
wmega = (exp(l2)*M[0,1]+exp(l1)*M[1,0])/(exp(l2)-exp(l1))

Wa = as_tensor([\
[0.0, wmega],\
[-wmega, 0.0 ]\
])

W = dot(dot(R,Wa),RT)

Wxy = conditional(gt(abs(l1-l2), 1.e-9), W[0,1], 0.0)
Wyx = conditional(gt(abs(l1-l2), 1.e-9), W[1,0], 0.0)

return (Bxx, Byx, Byy, Wxy, Wyx)

#---------------------------------------------------------------------------
# Function log_conformation
#---------------------------------------------------------------------------

def log_conformation( Sxx, Syx, Syy ):

# Calculate eigenvalues & eigenvectors of S tensor
(l1, l2, v1x, v1y) = eigenspectum( Sxx, Syx, Syy )

# Define R & R transpose tensors
R = as_tensor([\
[v1x, -v1y],\
[v1y, v1x]\
])

RT = transpose(R)

# Define conformation tensor
D = as_tensor([\
[exp(l1), 0.0 ],\
[0.0, exp(l2)]\
])

C = dot(dot(R,D),RT)

Cxx = C[0,0]
Cyx = C[0,1]
Cyy = C[1,1]

# Define e^(-S) tensor
Di = as_tensor([\
[exp(-l1), 0.0 ],\
[0.0, exp(-l2)]\
])

exp_mS = dot(dot(R,Di),RT)

exp_mSxx = exp_mS[0,0]
exp_mSyx = exp_mS[0,1]
exp_mSyy = exp_mS[1,1]

return (Cxx, Cyx, Cyy, exp_mSxx, exp_mSyx, exp_mSyy)

Community: FEniCS Project