Hello everyone, I bring you a question that has caused me great efforts without results. I need to apply a unidirectional stretch, however, as it is a vector element, which has two components, I do not know how to correctly deal with the matrices to be passed in emf. dirichletbc().
Although the code is relatively extensive, it presents an easy interpretation, however, due to the existence of EDLs (dielectric layers) and debye layer, it requires a high degree of refinement, therefore, I chose to use a mesh built and imported from gmsh , with a certain level of refinement and well-defined physical entities, instead of a basic MWE geometry.
Below I make the mesh and the code used available for those who are interested not only in helping me, but also for those who want some implementation reference for Boundary conditions, initial conditions, application of properties in subdomains, and the newmark discretization method.
mesh - Google Drive
Below is the stretching contour of the problem 2D.
It’s worth highlighting what I did for a two-dimensional case
Below is an implementation model that I tried to satisfy this condition.
from mpi4py import MPI
from dolfinx import fem,nls, log
import numpy as np
from ufl import VectorElement,FiniteElement,MixedElement,TestFunction,TrialFunction,split,grad,tr,Identity,inner,dot
from petsc4py.PETSc import ScalarType
import ufl
from dolfinx.io import gmshio
from petsc4py import PETSc
"""objective of the simulation"""
#study = "Electrochemical only"
study = "Capacitive strain sensing"
if study == "Electrochemical only":
T2_tot = 30.0e6
elif study == "Capacitive strain sensing":
T2_tot = 60.0e6
domain, cell_tags, facet_tags = gmshio.read_from_msh("malhas/capacitor_geometria__triangulo_13.msh", MPI.COMM_WORLD,0, gdim=2)
"""Function Space"""
# Define function space, scalar
U2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 1) # Displacent
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1) # concentrações + - , electric potential
D0 = ufl.FiniteElement("DG", domain.ufl_cell(), 0)
D1 = ufl.FiniteElement("DG", domain.ufl_cell(), 1)
# DOFs
TH = ufl.MixedElement([U2, P1, P1, P1])
ME = fem.FunctionSpace(domain, TH) # Total space for all DOFs
W = fem.FunctionSpace(domain,P1) # Scalar space for visualization later
"""Extraindo os sub espaços do elemento misto e os mapas contendo os graus de liberdade """
num_subs = ME.num_sub_spaces
spaces = []
maps = []
for i in range(num_subs):
space_i, map_i = ME.sub(i).collapse()
spaces.append(space_i)
maps.append(map_i)
"""Spaces for aplications """
"""
1 22 "wear_left"
1 23 "wear_right"
1 27 "ground"
1 28 "volt"
2 24 "Hidrogel_high"
2 25 "Hidrogel_low"
2 26 "elastomero"
"""
"""Extrair as tags da malha"""
def tag(n_tag):
return cell_tags.find(n_tag)
"Função para determinação das propriedades de um material em cada subdominio "
def mat_features(function_descontinuo, material, constanste):
space_Mat = fem.Function(function_descontinuo)
for i in range(len(material)):
space_Mat.x.array[material[i]] = np.full_like(material[i], constanste[i], dtype=ScalarType)
return space_Mat
'''''''''''''''''''''''''''''''''''''''''
Implementação das propriedades em cada regiao
'''''''''''''''''''''''''''''''''''''''
Q = fem.FunctionSpace(domain, D0)
hidrogel_sup = tag(24)
hidrogel_inf = tag(25)
elastomero = tag(26)
# A ordem de entrada das propriedades na lista deve ser equivalente ao espaço no qual o material ocupa dentro do dominio
material_ = [elastomero, hidrogel_inf, hidrogel_sup]
Gshear = mat_features(Q, material_, [0.034e-6, 0.003e-6, 0.003e-6])
Kbulk = mat_features(Q, material_, [2000*0.034e-6,2000*0.003e-6, 2000*0.003e-6])
#intInd = mat_features(Q, material_, [1, 0, 0])
matInd = mat_features(Q, material_, [0.0, 1, 1])
Gshear0 = 100.0e-6 # uso na formução fraca e no cconstrutor
Im_gent = mat_features(Q, material_, [90, 300.0, 300.0])
"""Constantes"""
D = 1.0e-2 # 1.0e0 # Diffusivity [L2T-1]
RT = 8.3145e-3 * (273.0 + 20.0) # Gas constant*Temp [ML2T-2#-1]
Farad = 96485.0e-6 # Faraday constant [Q#-1]
""" Initial concentrations"""
cPos0 = 0.274 # Initial concentration [#L-3]
cNeg0 = cPos0 # Initial concentration [#L-3]
cMax = 10000 * 1e-9 # 0.00001*1e-9 # 10000e0*1.e-9
"""Eletrical permittivity """
vareps0 = fem.Constant(domain, 8.85e-12 * 1e-6) #
vareps_num = mat_features(Q, material_, [1.0 ,1.0e3, 1.0e3]) # Permissividade do gel igual a 10000
vareps_r = mat_features(Q, material_, [6.5, 80., 80.]) # eletrical permittivity for material
vareps = vareps0 * vareps_r * vareps_num
# Mass density
rho = fem.Constant(domain, 1e-9) # 1e3 kg/m^3 = 1e-9 mg/um^3,
# Rayleigh damping coefficients
eta_m = 0.00000 # Constant(0.000005)
eta_k = 0.00000 # Constant(0.000005)
""" Generalized-alpha method parameters """
alpha = 0.0
gamma = 0.5 + alpha
beta = ((gamma + 0.5) ** 2) / 4
""" Simulation time related params (reminder: in microseconds)"""
ttd = 0.01
# Step in time
t = 0.0 # initialization of time
T_tot = 0e6 * 1.0e6 # 200.0 #0.11 # total simulation time
dt = T2_tot / 20 # incrementos de tempo
phi_norm = RT / Farad # "Thermal" Volt
"""unknown functions, testing functions and step functions"""
(u_test, omgPos_test, phi_test, omgNeg_test) = ufl.TestFunctions(ME) # Test function
# Define test functions in weak form
dw = ufl.TrialFunction(ME)
# Define actual functions with the required DOFs
w = fem.Function(ME)
(u, omgPos, phi, omgNeg) = ufl.split(w) # current DOFs
# A copy of functions to store values in last step for time-stepping.
w_old = fem.Function(ME)
(u_old, omgPos_old, phi_old, omgNeg_old) = ufl.split(w_old) # old DOFs
"""acceleration and speed"""
W2 = fem.FunctionSpace(domain, U2)
v_old = fem.Function(W2)
a_old = fem.Function(W2)
"""Implementação das Concentrações iniciais, cations e anions"""
x = ufl.SpatialCoordinate(domain)
DOlfin_Eps= 0.3e-15
# Initial electro-chemical potentials
Init_CPos= fem.Function(Q)
Init_CPos.x.array[material_[0]]= ufl.ln(DOlfin_Eps)
Init_CPos.x.array[material_[1]]= ufl.ln(cPos0)
Init_CPos.x.array[material_[2]]= ufl.ln(cPos0)
Init_CNeg= fem.Function(Q)
Init_CNeg.x.array[material_[0]]= ufl.ln(DOlfin_Eps)
Init_CNeg.x.array[material_[1]]= ufl.ln(cNeg0)
Init_CNeg.x.array[material_[2]]= ufl.ln(cNeg0)
w_old.sub(1).interpolate(Init_CPos)
w_old.sub(3).interpolate(Init_CNeg)
w.sub(3).interpolate(Init_CNeg)
w.sub(1).interpolate(Init_CPos)
'''''''''''''''''''''''''''''''''''''''''
BOUNDARY CONDITIONS
'''''''''''''''''''''''''''''''''''''''
fdim = domain.topology.dim - 1
def fixed_displacement_expression(x):
return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1])))
fixed_displacement = fem.Function(spaces[0])
fixed_displacement.interpolate(fixed_displacement_expression)
"""Engaste"""
Engast_0 = fem.locate_dofs_topological((ME.sub(0), spaces[0]), fdim, facet_tags.find(22))
Engast0 = fem.dirichletbc(fixed_displacement, Engast_0, ME.sub(0))
"""Estiramento : Deve-se alterar para obter os valores ideais de estiramento"""
scaleX = 2.0e4
class disp_Exp:
def __init__(self):
self.t= 0.0
self.Tramp =T2_tot
self.scalex = scaleX
def eval(self,x):
return np.stack(( np.full(x.shape[1], 5.5*self.scalex*self.t/self.Tramp),np.zeros(x.shape[1])))
def eval2(self,x):
return np.full(x.shape[1], 5.5*self.scalex*self.t/self.Tramp)
dispV= disp_Exp()
dispV.t= 20e6
dispV.Tramp= T2_tot
dispV.scalex= scaleX
disp= fem.Function(spaces[0])
disp.interpolate(dispV.eval)
#stretch00= fem.locate_dofs_topological((ME.sub(0),spaces[0]),fdim,facet_tags.find(23))
U_xy, Umaps= ME.sub(0).collapse()
Ux, sububsUx= U_xy.sub(0).collapse()
dispa= fem.Function(Ux)
dispa.interpolate(dispV.eval2)
stretch00= fem.locate_dofs_topological((ME.sub(0).sub(0),Ux),fdim,facet_tags.find(23))
bc_disp= fem.dirichletbc(dispa, stretch00, ME.sub(0))
#bc=[Engast0,bc_disp]
I made two attempts:
a) using the “conventional” way, which I found in some reference codes, where simply extract the X component through V.sub(0).sub(0).
b) creating a function with only the components of x, interpolating the function and applying it to fem.dirichletbc().
The problem is that the matrix I provide to
fem.dirichletbc()
It’s not suitable, because I make the Y component embedded, which I don’t want, however, if I put another argument, I get segmentation fault or matrix size faults. NOTE: The mesh is correct, given that the results of electrical potential and concentrations are in accordance with the literature, therefore, the problem is in the implementation of the stretching.
below I present the complete code:
from mpi4py import MPI
from dolfinx import fem,nls, log
import numpy as np
from ufl import VectorElement,FiniteElement,MixedElement,TestFunction,TrialFunction,split,grad,tr,Identity,inner,dot
from petsc4py.PETSc import ScalarType
import ufl
from dolfinx.io import gmshio
from petsc4py import PETSc
domain, cell_tags, facet_tags = gmshio.read_from_msh("malhas/capacitor_geometria__triangulo_13.msh", MPI.COMM_WORLD,0, gdim=2)
"""Function Space"""
# Define function space, scalar
U2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 1) # Displacent
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1) # concentrações + - , electric potential
D0 = ufl.FiniteElement("DG", domain.ufl_cell(), 0)
D1 = ufl.FiniteElement("DG", domain.ufl_cell(), 1)
# DOFs
TH = ufl.MixedElement([U2, P1, P1, P1])
ME = fem.FunctionSpace(domain, TH) # Total space for all DOFs
W = fem.FunctionSpace(domain,P1) # Scalar space for visualization later
"""Extraindo os sub espaços do elemento misto e os mapas contendo os graus de liberdade """
num_subs = ME.num_sub_spaces
spaces = []
maps = []
for i in range(num_subs):
space_i, map_i = ME.sub(i).collapse()
spaces.append(space_i)
maps.append(map_i)
"""Aplicando as propriedades dos materiais"""
"""
1 28 "Voltagem"
1 27 "Aterramento"
1 22 "Engaste_esquerda"
1 23 "Engaste_direita"
2 24 "Hidrogel_superior"
2 25 "Hidrogel_inferior"
2 26 "elastomero"
2 27 "dominio"
"""
"""Extrair as tags da malha"""
def tag(n_tag):
return cell_tags.find(n_tag)
"Determinação das propriedades de um material "
def mat_features(function_descontinuo, material, constanste):
space_Mat = fem.Function(function_descontinuo)
for i in range(len(material)):
space_Mat.x.array[material[i]] = np.full_like(material[i], constanste[i], dtype=ScalarType)
return space_Mat
""" Implementação das propriedades em cada regiao"""
Q = fem.FunctionSpace(domain, D0)
hidrogel_sup = tag(24)
hidrogel_inf = tag(25)
elastomero = tag(26)
# A ordem de entrada das propriedades na lista deve ser equivalente ao espaço no qual o material ocupa dentro do dominio
material_ = [elastomero, hidrogel_inf, hidrogel_sup]
Gshear = mat_features(Q, material_, [0.034e-6, 0.003e-6, 0.003e-6])
Kbulk = mat_features(Q, material_, [2000*0.034e-6,2000*0.003e-6, 2000*0.003e-6])
#intInd = mat_features(Q, material_, [1, 0, 0])
matInd = mat_features(Q, material_, [0.0, 1, 1])
Gshear0 = 100.0e-6 # uso na formução fraca e no cconstrutor
Im_gent = mat_features(Q, material_, [90, 300.0, 300.0])
"""Constantes"""
D = 1.0e-2 # 1.0e0 # Diffusivity [L2T-1]
RT = 8.3145e-3 * (273.0 + 20.0) # Gas constant*Temp [ML2T-2#-1]
Farad = 96485.0e-6 # Faraday constant [Q#-1]
""" Initial concentrations"""
cPos0 = 0.274 # Initial concentration [#L-3]
cNeg0 = cPos0 # Initial concentration [#L-3]
cMax = 10000 * 1e-9 # 0.00001*1e-9 # 10000e0*1.e-9
"""Eletrical permittivity """
vareps0 = fem.Constant(domain, 8.85e-12 * 1e-6) #
vareps_num = mat_features(Q, material_, [1.0 ,1.0e3, 1.0e3]) # Permissividade do gel igual a 10000
vareps_r = mat_features(Q, material_, [6.5, 80., 80.]) # eletrical permittivity for material
vareps = vareps0 * vareps_r * vareps_num
# Mass density
rho = fem.Constant(domain, 1e-9) # 1e3 kg/m^3 = 1e-9 mg/um^3,
# Rayleigh damping coefficients
eta_m = 0.00000 # Constant(0.000005)
eta_k = 0.00000 # Constant(0.000005)
""" Generalized-alpha method parameters """
alpha = 0.0
gamma = 0.5 + alpha
beta = ((gamma + 0.5) ** 2) / 4
""" Simulation time related params (reminder: in microseconds)"""
ttd = 0.01
# Step in time
t = 0.0 # initialization of time
T_tot = 0e6 * 1.0e6 # 200.0 #0.11 # total simulation time
study = "Electrochemical only"
#study = "Capacitive strain sensing"
if study == "Electrochemical only":
T2_tot = 30.0e6
elif study == "Capacitive strain sensing":
T2_tot = 60.0e6
dt = T2_tot / 20 # incrementos de tempo
phi_norm = RT / Farad # "Thermal" Volt
"""unknown functions, testing functions and step functions"""
(u_test, omgPos_test, phi_test, omgNeg_test) = ufl.TestFunctions(ME) # Test function
# Define test functions in weak form
dw = ufl.TrialFunction(ME)
# Define actual functions with the required DOFs
w = fem.Function(ME)
(u, omgPos, phi, omgNeg) = ufl.split(w) # current DOFs
# A copy of functions to store values in last step for time-stepping.
w_old = fem.Function(ME)
(u_old, omgPos_old, phi_old, omgNeg_old) = ufl.split(w_old) # old DOFs
"""acceleration and speed"""
W2 = fem.FunctionSpace(domain, U2)
v_old = fem.Function(W2)
a_old = fem.Function(W2)
"""Implementação das Concentrações iniciais, cations e anions, initial concentrations"""
x = ufl.SpatialCoordinate(domain)
DOlfin_Eps= 0.3e-15
# Initial electro-chemical potentials
Init_CPos= fem.Function(Q)
Init_CPos.x.array[material_[0]]= ufl.ln(DOlfin_Eps)
Init_CPos.x.array[material_[1]]= ufl.ln(cPos0)
Init_CPos.x.array[material_[2]]= ufl.ln(cPos0)
Init_CNeg= fem.Function(Q)
Init_CNeg.x.array[material_[0]]= ufl.ln(DOlfin_Eps)
Init_CNeg.x.array[material_[1]]= ufl.ln(cNeg0)
Init_CNeg.x.array[material_[2]]= ufl.ln(cNeg0)
w_old.sub(1).interpolate(Init_CPos)
w_old.sub(3).interpolate(Init_CNeg)
w.sub(3).interpolate(Init_CNeg)
w.sub(1).interpolate(Init_CPos)
'''''''''''''''''''''''''''''''''''''''''
BOUNDARY CONDITIONS
'''''''''''''''''''''''''''''''''''''''
fdim = domain.topology.dim - 1
"""Engaste (wear) """
def fixed_displacement_expression(x):
return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1])))
fixed_displacement = fem.Function(spaces[0])
fixed_displacement.interpolate(fixed_displacement_expression)
Engast_0 = fem.locate_dofs_topological((ME.sub(0), spaces[0]), fdim, facet_tags.find(22))
Engast0 = fem.dirichletbc(fixed_displacement, Engast_0, ME.sub(0).sub(0))
"""Estiramento : Deve-se alterar para obter os valores ideais de estiramento"""
scaleX = 2.0e4
class disp_Exp:
def __init__(self):
self.t= 0.0
self.Tramp =T2_tot
self.scalex = scaleX
def eval(self,x):
return np.stack(( np.full(x.shape[1], 5.5*self.scalex*self.t/self.Tramp),np.zeros(x.shape[1])))
dispV= disp_Exp()
dispV.t= 0.
dispV.Tramp= T2_tot
dispV.scalex= scaleX
disp= fem.Function(spaces[0])
disp.interpolate(dispV.eval)
stretch00= fem.locate_dofs_topological((ME.sub(0),spaces[0]),fdim,facet_tags.find(23))
bc_disp= fem.dirichletbc(disp, stretch00, ME.sub(0).sub(0))
"""Aterramento"""
def ground_0(x):
return np.stack((np.zeros(x.shape[1])))
ground = fem.Function(spaces[2])
ground.interpolate(ground_0)
ground0 = fem.locate_dofs_topological((ME.sub(2), spaces[2]), fdim, facet_tags.find(27))
bc_ground = fem.dirichletbc(ground, ground0, ME.sub(2))
"""Aplicação da voltagem inicial de 0 à 1, em um intervalo de 30s, e permanecer constante em 1 indefinidamente """
class phiRamp_function():
def __init__(self):
self.phi_norm= phi_norm
self.pi=np.pi
self.Tramp= 30.0e6
self.t= 0.0
self.temp= 1.0e3
def phi_eval(self,x):
return np.stack((np.full(x.shape[1], min((self.temp/self.phi_norm)*(self.t/self.Tramp), self.temp/self.phi_norm))))
phiRamp_func= phiRamp_function()
phiRamp_func.t= 0
phiRamp= fem.Function(spaces[2])
phiRamp.interpolate(phiRamp_func.phi_eval)
phiRamp_0= fem.locate_dofs_topological((ME.sub(2),spaces[2]),fdim,facet_tags.find(28))
bc_phiRamp = fem.dirichletbc(phiRamp,phiRamp_0,ME.sub(2))
bc=[Engast0, bc_ground,bc_phiRamp,bc_disp]
'''''''''''''''''''''
SUBROUTINES
'''''''''''''''''''''
dk = fem.Constant(domain,dt)
def F_calc(u):
dim = len(u)
Id = Identity(dim) # Identity tensor
F = Id + grad(u) # 2d Deformation gradient
return F # Full 2d F
def Piola(F,phi):
Id = Identity(2)
J = ufl.det(F)
C = F.T*F
Cdis = J**(-2/3)*C
I1 = tr(Cdis)
eR = -phi_norm*grad(phi)
e_sp = ufl.inv(F.T)*eR
#
T_max = vareps0*vareps_r*J*(ufl.outer(e_sp,e_sp) - 1/2*(inner(e_sp,e_sp))*Id)*ufl.inv(F.T)
# Piola stress (Gent)
TR = J**(-2/3)*Gshear*(Im_gent/(Im_gent+3-I1))*(F - 1/3*tr(Cdis)*ufl.inv(F.T))\
+ Kbulk*ufl.ln(J)*ufl.inv(F.T) + T_max
return TR
# Update formula for acceleration
# a = 1/(2*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) - (1-2*beta)*a0)
def update_a(u, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dk
beta_ = beta
else:
dt_ = float(dk)
beta_ = float(beta)
return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old
# Update formula for velocity
# v = dt * ((1-gamma)*a0 + gamma*a) + v0
def update_v(a, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dk
gamma_ = gamma
else:
dt_ = float(dk)
gamma_ = float(gamma)
return v_old + dt_*((1-gamma_)*a_old + gamma_*a)
def update_fields(u_proj, u_proj_old, v_old, a_old):
"""Update fields at the end of each time step."""
# Get vectors (references)
u_vec, u0_vec = u_proj.x.array[:], u_proj_old.x.array[:]
v0_vec, a0_vec = v_old.x.array[:], a_old.x.array[:]
# use update functions using vector arguments
a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl=False)
v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl=False)
# Update (u_old <- u)
v_old.x.array[:], a_old.x.array[:] = v_vec, a_vec
#v_old.vector()[:], a_old.vector()[:] = v_vec, a_vec
v_old.x.scatter_forward()
a_old.x.scatter_forward()
def avg(x_old, x_new, alpha):
return alpha*x_old + (1-alpha)*x_new
'''''''''''''''''''''''''''''''''''''''''
KINEMATICS & CONSTITUTIVE RELATIONS
'''''''''''''''''''''''''''''''''''''''''
# Newmark-beta kinematical update
a_new = update_a(u, u_old, v_old, a_old, ufl=True)
v_new = update_v(a_new, u_old, v_old, a_old, ufl=True)
# get avg fields for generalized-alpha method
u_avg = avg(u_old, u, alpha)
v_avg = avg(v_old, v_new, alpha)
omgPos_avg = avg(omgPos_old, omgPos, alpha)
phi_avg = avg(phi_old, phi, alpha)
omgNeg_avg = avg(omgNeg_old, omgNeg, alpha)
# Explicit concentration updates
cPos = ufl.exp(omgPos_avg - phi_avg)
cNeg = ufl.exp(omgNeg_avg + phi_avg)
cPos_old = ufl.exp(omgPos_old - phi_old)
cNeg_old = ufl.exp(omgNeg_old + phi_old)
# Kinematics
F = F_calc(u_avg)
C = F.T*F
Ci = ufl.inv(C)
F_old = F_calc(u_old)
J = ufl.det(F)
J_old = ufl.det(F_old)
'''''''''''''''''''''''
WEAK FORMS
'''''''''''''''''''''''
# Enable inertial forces
dynSwitch = fem.Constant(domain, 1.0)
#dx = ufl.Measure('dx', domain=mesh, subdomain_data=materials)
dx = ufl.Measure('dx', domain=domain)
# Weak forms
L0 = inner(Piola(F_calc(u_avg), phi_avg), grad(u_test))*dx \
+ dynSwitch*rho*inner(a_new, u_test)*dx
L1 = dot((cPos-cPos_old)/dk/D,omgPos_test)*dx \
+ (cPos_old)*inner(Ci*matInd*grad(omgPos_avg),grad(omgPos_test))*dx
L2 = dot(vareps*phi_norm*J*Ci*grad(phi_avg),grad(phi_test))*dx \
-dot(Farad*matInd*cMax*(cPos-cNeg),phi_test)*dx
L3 = dot((cNeg-cNeg_old)/dk/D,omgNeg_test)*dx \
+ (cNeg_old)*dot(Ci*matInd*grad(omgNeg_avg),grad(omgNeg_test))*dx
# Total weak form
L = (1/Gshear0)*L0 + L1 + L2 + L3
# Automatic differentiation tangent:
a = ufl.derivative(L, w, dw)
problem = fem.petsc.NonlinearProblem(L,w,bc,J=a)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
log.set_log_level(log.LogLevel.INFO)
solver.atol = 1e-3
solver.rtol = 1e-3
solver.convergence_criterion = "incremental"
solver.max_it =1000
solver.report = True
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
from dolfinx.io import XDMFFile
from datetime import datetime
file0 = XDMFFile(MPI.COMM_WORLD, "resultados/capacitor0.xdmf", "w")
file0.write_mesh(domain)
file1 = XDMFFile(MPI.COMM_WORLD, "resultados/capacitor1.xdmf", "w")
file1.write_mesh(domain)
file2 = XDMFFile(MPI.COMM_WORLD, "resultados/capacitor2.xdmf", "w")
file2.write_mesh(domain)
file3 = XDMFFile(MPI.COMM_WORLD, "resultados/capacitor3.xdmf", "w")
file3.write_mesh(domain)
file4 = XDMFFile(MPI.COMM_WORLD, "resultados/capacitor4.xdmf", "w")
file4.write_mesh(domain)
"""saves"""
phi_volt= fem.Function(W)
def writeResults(t):
# Convert phi to Volts before saving
phi_Dim = phi_norm/1.e3*phi
phi_f= fem.Expression(phi_Dim,W.element.interpolation_points())
phi_volt.interpolate(phi_f)
file0.write_function(w_old.sub(0),t)
file1.write_function(w_old.sub(1),t)
file2.write_function(w_old.sub(2),t)
file3.write_function(w_old.sub(3),t)
file4.write_function(phi_volt,t)
"""Solver"""
while (round(t,2) <= round(T2_tot,2)):
writeResults(t)
# Output storage, also outputs intial state at t=0
if t<=30.0e6:
dispV.t = 0.0
else:
dispV.t = t - (alpha*dt) - 30.0e6
phiRamp_func.t = t - (alpha*dt)
phiRamp.interpolate(phiRamp_func.phi_eval)
disp.interpolate(dispV.eval)
# Solve the problem
iter, converged = solver.solve(w)
# Update fields for next step
u_proj=fem.Function(W2)
u_proj1 = fem.Expression(u, W2.element.interpolation_points())
u_proj.interpolate(u_proj1)
u_proj_old=fem.Function(W2)
u_proj_old1 =fem.Expression(u_old, W2.element.interpolation_points())
u_proj_old.interpolate(u_proj_old1)
update_fields(u_proj, u_proj_old, v_old, a_old)
w_old.x.array[:] = w.x.array
now = datetime.now()
current_time = now.strftime("%H:%M:%S")
print("Step: Sense | Simulation Time: {} s | Iterations: {}".format(t/1e6, iter))
print()
# increment time, counters
t += dt
print(t)
#last step
file0.write_function(w_old.sub(0),t)
file1.write_function(w_old.sub(1),t)
file2.write_function(w_old.sub(2),t)
file3.write_function(w_old.sub(3),t)```