Initial Conditions and Function update on time in mixed Element

Version: Dolfinx 0.6.0 Python
Windows Subsystem for Linux (WSL) in Visual Studio Code

Hello guys, I’ve been trying to solve an Eletro-chemo-mechanical coupling problem, and I’ve been faced with some problems.
I’m trying to update these codes to the current version FENICSX, however, I’m having a lot of difficulties.
Due to the fact that the code is extremely extensive, I chose not to write it in full here, just to present the fundamental points

Some questions are:

  1. How to determine an initial concentration the of mobile ions c0= 2740 , in time t =0 for each subdomain ( hydrogel and elastomer) , for a mixed elements.

  2. How to update scalar and vector solution functions correctly in newmark functions (displacement, speed and acceleration)

In the previous version,
This was done in the following way

from dolfin import *
import numpy as np

# Define function space, scalar
U2 = VectorElement("Lagrange", mesh.ufl_cell(), 1)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
D0 = FiniteElement("DG", mesh.ufl_cell(), 0)
D1 = FiniteElement("DG", mesh.ufl_cell(), 1)

# DOFs
TH = MixedElement([U2, P1, P1, P1])
ME = FunctionSpace(mesh, TH) # Total space for all DOFs

W = FunctionSpace(mesh,P1)   # Scalar space for visualization later
W2 = FunctionSpace(mesh,U2)   # Vector space for visualization later
W3 = FunctionSpace(mesh,D0)   # DG space for visualization later
W4 = FunctionSpace(mesh,D1)   # DG space for visualization later

# Define test functions in weak form
dw = TrialFunction(ME)                                   
(u_test, omgPos_test, phi_test, omgNeg_test)  = TestFunctions(ME)    # Test function

# Define actual functions with the required DOFs
w = Function(ME)
(u, omgPos, phi, omgNeg) = split(w)    # current DOFs

# A copy of functions to store values in last step for time-stepping.
w_old = Function(ME)
(u_old, omgPos_old, phi_old, omgNeg_old) = split(w_old)   # old DOFs

v_old = Function(W2)
a_old = Function(W2)

# Initial electro-chemical potentials
init_omgPos = Expression('abs(x[1])>=int2-tol?std::log((cPos0)):std::log((cNum))', int2=int2, tol = tol, cPos0 = cPos0, cNum=DOLFIN_EPS, degree=0)

omgPos_init = interpolate(init_omgPos,ME.sub(1).collapse())
assign(w_old.sub(1),omgPos_init)
assign(w.sub(1),omgPos_init)


# 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.vector(), u_proj_old.vector()
    v0_vec, a0_vec = v_old.vector(), a_old.vector()

    # 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.vector()[:], a_old.vector()[:] = v_vec, a_vec
    #u_old.vector()[:] = u_proj.vector()

def avg(x_old, x_new, alpha):
    return alpha*x_old + (1-alpha)*x_new


# Update fields for next step
    u_proj = project(u, W2)
    u_proj_old = project(u_old, W2)
    update_fields(u_proj, u_proj_old, v_old, a_old)
    w_old.vector()[:] = w.vector()
    

I made some attempts to solve the problems, however I am not making any progress in terms of updating the vector fields and applying the initial condition, due to the fact that I must determine the initial concentration in different subdominions, for a geometry built in gmsh.
I made two attempts to determine the concentrations, one from the coordinates and the other through physical entities built in gmsh itself.Below is the summarized code.

From the coordinates:

from mpi4py import MPI
from dolfinx import fem, io, nls, log, mesh, plot
import numpy as np
import pyvista
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

domain, cell_tags, facet_tags = gmshio.read_from_msh("capacitor_malha.msh", MPI.COMM_SELF,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

"""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)



(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

"""Concentrações iniciais"""
x = ufl.SpatialCoordinate(domain)

DOlfin_Eps= 0.3e-16
scaleY = 900.e0 #Problema aqui
tol = 1e-12
cPos0 = 0.274
int1 = 200.e0 # 
int2 = scaleY/2.-int1

class init_Omg:
    def __init__(self):
        self.scaley= scaleY
        self.int2= int2
        self.tol= tol
        self.c= 0
        self.CNum= DOlfin_Eps
        
    def init_Omg_eval(self,x):
        if x[1]>= 700e-6 + tol:
            #return ufl.ln(self.c) #
            return np.log(self.c, dtype=np.int32)
        if x[1]<= 200e-6 - tol:
            return np.log(self.c, dtype=np.int32)
        else:
            return np.log(self.CNum, dtype=np.int32)

    def init_Omg_(self,x):
        return ufl.conditional(ufl.ge(x[1],200e+6 + tol), ufl.ln(self.c), ufl.ln(self.CNum))
    


my_bc_value = init_Omg()
my_bc_value.c=cPos0
bc_expr = fem.Expression(my_bc_value.init_Omg_eval, spaces[1].element.interpolation_points())
w_old.sub(1).interpolate(bc_expr)

by physical entities:


""" Implementation of concentrations in each region"""
def tag(n_tag):
    return cell_tags.find(n_tag)

Q = fem.FunctionSpace(domain, D0)

eletrodo_sup = tag(26)
eletrodo_inf = tag(27)
gel_p = tag(28)

# A ordem de entrada das propriedades na lista deve ser equivalente ao espaço no qual o material ocupa dentro do dominio
material_ = [gel_p, eletrodo_inf, eletrodo_sup]

Init_CPos= fem.Function(spaces[1])
Init_CPos.interpolate(my_bc_value.init_Omg_)

w_old.x.array[gel_p]= Init_CPos.x.array

both attempts don’t work.

I guess you are referring to this?

In the code below, you are mixing a whole lot of ufl and numpy

Let us consider a simple MWE, on a unit square:

from IPython import embed
from mpi4py import MPI
from dolfinx import fem, io, nls, log, mesh, plot
import numpy as np
from petsc4py.PETSc import ScalarType
import ufl
from dolfinx.io import gmshio


domain = mesh.create_rectangle(MPI.COMM_WORLD, [[2, 2], [3, 3]], [10, 10])
"""Function Space"""
# Define function space, scalar
U2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 1)  # Displacent
# concentrações + - , electric potential
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
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.Function(ME)

a = 2.5
cpos0 = 1.1
cNum = 2


def init_omega(x):
    condition = x[1] > a
    condition_inverse = np.invert(condition)
    return condition*np.log(cpos0) + condition_inverse * np.log(cNum)
w.sub(1).interpolate(init_omega)

Here I’ve created a simply numpy function with the same functionality as

where you replace a with int2-tol etc.

As long as your example is kept as partial snippets, it is quite hard to give you any further guidance (and it is aso the reason why it has taken so long to write up a reply). I would strongly suggest you try to simplify your code, and have a look at: PETSc NonlinearProblem running slow - #4 by dokken
as it seems rather relevant for the second question.