Source term as an Expression

variational formulation
1

Hello everybody,

i have the following problem.
I want to implement my source term f as an Expression.

gamma_theta = 0.2. # electric conductivity
eta = 0.24                 # heat conductivity

f = Expression("eta*pow(abs(S))+gamma_theta*pow(abs(grad(elektrisches_Pot)))", eta=eta,S=s,gamma_theta=gamma_theta,elektrisches_Pot=elektrisches_Pot)

S and elektrisches_Pot are calculated with two previous problems like

solve(a == L, u, bc2)
s = sigma(u_h2) - (1./2)*tr(sigma(u_h2))*Identity(d)  # deviatoric stress in 2D

where u is the solution of the equation of potential.

I am getting the following error:

RuntimeError: Unable to compile C++ code with dijitso

It would be really nice if someone could help me out here.
Thank you!

2

Here is a MWE

from dolfin import *
import numpy as np

mesh = UnitSquareMesh(10, 10)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()

# Initialize mesh function for boundary domains
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)

boundary_subdomains.set_all(4)
left.mark(boundary_subdomains, 1)
top.mark(boundary_subdomains, 2)
right.mark(boundary_subdomains, 3)
bottom.mark(boundary_subdomains, 4)

degree = 1

dss=Measure('ds', domain=mesh, subdomain_data=boundary_subdomains)
area_gamma1 = assemble(1*dss(4))
# specific electric conductivity of carbon steel at 20 C
leitfaehigkeit = 1.43e-7

V_phi = FunctionSpace(mesh, "CG", degree)

# Dirichlet boundary condition for the equation of potential
c = Constant(0.0)
bc2 = DirichletBC(V_phi, c, boundary_subdomains, 2)
dx = Measure('dx', mesh)

# Test and trial functions for equation of potential
phi = TrialFunction(V_phi)
v = TestFunction(V_phi)

u_h1 = Function(V_phi)  # for equation of potential

f = Constant(0.0)

# bilinearform for the equation of potentialĂź
a = dot(grad(phi), grad(v))*dx


# time parameters
T = 2.0            # final time
num_steps = 10     # number of time steps
dt = T / num_steps # time step size
t = 0.0


alpha = 3          # parameter alpha  # those parameters are for now for the L2-projection and need to be changed
beta = 1.2         # parameter beta
        
# repleasment functions for the equation of potential in the discrete space
u = Function(V_phi)
u_n = Function(V_phi)

# second problem -------------------------------------------------------------------------------------------
V_stress = VectorFunctionSpace(mesh, "CG", degree)

# Test and trial functions for second problem
u2 = TrialFunction(V_stress)
v2 = TestFunction(V_stress)

# Lame coefficients
# 2D case for the creep strain-----------------------
E = Constant(1e5)  # Youngs module
nu = Constant(0.3)  # 
model = "plane_stress"

mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
    lmbda = 2*mu*lmbda/(lmbda+2*mu)


def eps(v):
    return sym(grad(v))

# stress tensor
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

faktor = Constant(0.5) # in 2D ist der Vorfaktor 1/n = 0.5

f2 = Constant((0, 0))  # no pressure were supplied
#T = Constant((1, 1))
T = np.random.rand(num_steps, 2)
print("Random stress on the boundary = " + str(T))
a2 = inner(sigma(u2), eps(v2))*dx

# Needed to avoid error in the linear form (see ufl datatypes)
T_ufl = as_matrix(T.tolist())

c2 = Constant((0.0, 0.0))
bc2_2 = DirichletBC(V_stress, c2, boundary_subdomains, 2)

u_h2 = Function(V_stress)  # for second problem
u_h2n = Function(V_stress)
d = u_h2.geometric_dimension()  # space dimension


# Problem heat equation-----------------------------------------------
V_heat = VectorFunctionSpace(mesh, "CG", 1)

# Define test and trial function
u_heat = TrialFunction(V_heat)
v_heat = TestFunction(V_heat)

# define \gamma(\theta), \eta(\theta) 
gamma_theta = 0.2
eta = Constant(0.24) # for \eta(\theta) heat conductivity

c3 = Constant((0.0, 0.0))
bc2_3 = DirichletBC(V_heat, c2, boundary_subdomains, 2)

u_h3 = Function(V_heat)  # for second problem
u_h3n = Function(V_heat)

sol = Constant(0.2345)  # just for testing; it is the solution from another class


for i in range(num_steps-1):

    # potential of equation-------------------------
    L = f*v*dx + (leitfaehigkeit*(1/area_gamma1)*(-sol))*v*dss(4)
    # Compute solution
    solve(a == L, u, bc2)
    u_n.assign(u)


    # stress stuff ---------------------------------------------------
    L2 = dot(f2, v2)*dx + dot(T_ufl[i,:], v2)*dss(4)
    solve(a2 == L2, u_h2, bc2_2)
       

    # stress
    s = sigma(u_h2) - (1./2)*tr(sigma(u_h2))*Identity(d)  # deviatoric stress in 2D
    von_Mises = sqrt(3./2*inner(s, s))
    V = FunctionSpace(mesh, 'P', 1)
    von_Mises = project(von_Mises, V)



    # Compute magnitude of displacement
    u_magnitude = sqrt(dot(u_h2, u_h2))
    u_magnitude = project(u_magnitude, V)

    # non-stationary heat equation -----------------------------------------
    f = Expression("eta*pow(abs(S))+gamma_theta*pow(abs(grad(elektrisches_Pot)))", eta=eta,S=s,gamma_theta=gamma_theta,elektrisches_Pot=u)

    t += dt
3

Without looking at the MWE, doesn’t pow() require two positional arguments? You have not provided the equation you’re trying to replicate but I think you are only providing pow() the base, not the exponent.

2 Likes
4

Hello Sven,

thank you for your reply.

Do you mean something like this

    f_heat = Expression("eta*pow(abs(S),2)+gamma_theta*pow(abs(grad(elektrisches_Pot)),2)", degree=1, eta=0.24,S=s,gamma_theta=0.2,elektrisches_Pot=u)

and whereby the second argument of pow defines the exponent 2?
Unfortunately this still gives me the error.

Greetings

5

There were a couple of mistakes here. I fixed it now.
Again, thank you for your help Sven.

6

Ok. I haven´t fixed it :(.
It is open again. I would be really grateful for your help.
Maybe i have no other choice then defining f through a user defined expression. Wish it wouldn’t be like that ;).
But i have problems understanding the concept. Is there anywhere a good explanation please?

7

Hello everybody,

I had a little bit of time at my other work and i am doing it now with the User Expression. See please here

            class rhs_heat(UserExpression):
                def __init__(self, gamma_theta, eta, s, **kwargs):
                    super().__init__(**kwargs)
                    self.gamma_theta = gamma_theta
                    self.eta = eta
                    self.s = s

                #Here is the function where the solution u should be used  
                def eval(self, value, x):
                    p = eta*pow(abs(s),2)+gamma_theta*pow(abs(grad(u)),2)
                    pdb.set_trace()
                    value [0] =  p

            f_heat6 = rhs_heat(gamma_theta=gamma_theta,eta=eta,s=s)

But for sure i am getting the following error:

ValueError: setting an array element with a sequence.

Is it possible to use the user expression with the shapes of abs(grad(u)) and abs(s)?

Thank you!

8

Please make a minimal example that can reproduce the error, as illustrated by the links in Read before posting: How do I get my question answered? - #4 by nate.

Please also read through the guidelines in the first post of the link above.

9

As dokken is advising, it’s hard to provide guidance without seeing the whole picture. You also don’t need to include extra baggage such as pdb. You ask

In your most recent post you’ve removed the fact that s is a tensor and u is a scalar. It would appear from your previous posts that:

All of this information is useful to know and has to piecemealed together unless you provide a concise MWE. Then it is easier to be on the same page.

When you call eval() in your UserExpression, you are expecting to return a scalar for your class rhs_heat. However you are setting rhs_heat with

f_heat6 = rhs_heat(gamma_theta=gamma_theta,eta=eta,s=s)

If you print p, it likely won’t be a scalar. I believe at every point x you are returning p evaluated at every point in your mesh.

Also, you said earlier

Do you want feedback on the original approach or this current one? You did not provide the errors from the original post that led you to this new approach. Clarifying your post will help everyone get you an answer faster with less energy.

1 Like
10

Thank you Dokken and Sven for your useful comments.
I am sorry for not writing that s is a tensor and u is a scalar. It’s true, that would have been useful.

I solved it now (hopefully once and for all) with the value_shape.