# Temperature dependent material data from file

**URL:** <https://fenicsproject.discourse.group/t/temperature-dependent-material-data-from-file/4975>\
**Category:** Uncategorized\
**Created:** [January 26, 2021, 2:28pm UTC](https://fenicsproject.discourse.group/t/temperature-dependent-material-data-from-file/4975 "2021-01-26T14:28:29Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![mhechthz](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mhechthz/32/2100_2.png) [@mhechthz](https://fenicsproject.discourse.group/u/mhechthz)\
**Post date:** [January 26, 2021, 2:28pm UTC](https://fenicsproject.discourse.group/t/temperature-dependent-material-data-from-file/4975/1 "2021-01-26T14:28:29Z")

</div>

Hi,

for my time dependent heat convection-diffusion equation I have some material data from files: density, conductivity and specific heat. All depend on the current local temperature, i.e. local solution of the PDE. I tried the following, that works well with constant material data (as `Constant()` or `Expression()`):

```
...

# material data from file, int(temperature) == index
cpX = np.genfromtxt('analyse/cpX.ana')[:,1]
rhoX = np.genfromtxt('analyse/densityX.ana')[:,1]
lmbdaX = np.genfromtxt('analyse/conductX.ana')[:,1]

...

class userMaterial(UserExpression):

    def __init__ (self, temp, valArray, **kwargs):
        super(). __init__ (kwargs)
        self.temp = temp
        self.valArray = valArray

    def eval(self, value, x):
        intT = max(1,min(int(self.temp(x)),1600))
        value[0] = self.valArray[intT]

    def value_shape(self):
        return ()

...

U = Function(V)
U_data = np.array(U.vector())*0+100
U.vector().set_local(U_data)

lambdaU = userMaterial(temp=U,valArray=lmbdaX,degree=0)
rhoU = userMaterial(temp=U,valArray=rhoX,degree=0)
cpU = userMaterial(temp=U,valArray=cpX,degree=0)

...

E_thd = lambdaU/(rhoU*cpU) # thermal diffusivity
E_rcp = 1.0/(rhoU*cpU) # 1 / rho * cp

...

for i in bcond:
    if 'Dirichlet' in bcond[i]:
        bc = DirichletBC(V,bcond[i]['Dirichlet'],boundary_markers,i)
        bcs.append(bc)

    elif 'Neumann' in bcond[i]:
        if bcond[i]['Neumann'] != 0:
            g = bcond[i]['Neumann']
            integrals_N.append(E_dt*E_rcp*g*v*ds(i))

    elif 'Robin' in bcond[i]:
        r,s = bcond[i]['Robin']
        integrals_R_a.append(E_dt*E_rcp*r*u*v*ds(i))
        integrals_R_L.append(E_dt*E_rcp*r*s*v*ds(i))

...

#===========================================================
# PDE
#
# rho*cp*grad(T,t)+rho*cp*w*grad(T,x)=div(lambda*grad(T,x),x)
#
#===========================================================

# weak problem formulation

F = u*v*dx + E_dt*E_thd*dot(grad(u),grad(v))*dx - u_n*v*dx \
       + E_dt*inner(E_w,grad(u))*v*dx \
       + sum(integrals_R_a) + sum(integrals_N) - sum(integrals_R_L)

a, L = lhs(F), rhs(F)

...

t = 0
while t < tFinal:

    # Update current time
    t += dt
    print(t,"s /",tFinal,'s')

    # Compute solution
    solve(a == L, U, bcs)

    # Save to file and plot solution
    vtkfile << (U, t)
    # plot(u)

    # Update previous solution
    u_n.assign(U)

    # update material data
    lambdaU.temp = U
    rhoU.temp = U
    cpU.temp = U

...

```

Unfortunately, this doesn’t work and is also extremely slow.

How can measured, temperature dependend table data be applied to the PDE?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 26, 2021, 4:33pm UTC](https://fenicsproject.discourse.group/t/temperature-dependent-material-data-from-file/4975/2 "2021-01-26T16:33:58Z")

</div>

Write it into a function once, save it as a checkpoint using `XDMFFile.write_checkpoint` such that it can be reloaded in parallel with basically no cost using `XDMFFile.read_checkpoint` for later experiments.

---

<div class="post-metadata">

**Author:** ![mhechthz](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mhechthz/32/2100_2.png) [@mhechthz](https://fenicsproject.discourse.group/u/mhechthz)\
**Post date:** [January 27, 2021, 7:55am UTC](https://fenicsproject.discourse.group/t/temperature-dependent-material-data-from-file/4975/3 "2021-01-27T07:55:46Z")

</div>

Wow, this sounds as you know what you do. But I don’t understand much, since I’m new to FEniCS ☹

Could you add some pseudo code?

What I thought meanwhile is to use a MeshFunction and update it only on the vertices to avoid the expencive call to temp(x). Would this be possible?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 27, 2021, 7:57am UTC](https://fenicsproject.discourse.group/t/temperature-dependent-material-data-from-file/4975/4 "2021-01-27T07:57:03Z")

</div>

See for instance: [Read mesh from XDMF file (write\_checkpoint) - #3 by dokken](https://fenicsproject.discourse.group/t/read-mesh-from-xdmf-file-write-checkpoint/3458/3)

---

<div class="post-metadata">

**Author:** ![mhechthz](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mhechthz/32/2100_2.png) [@mhechthz](https://fenicsproject.discourse.group/u/mhechthz)\
**Post date:** [January 28, 2021, 2:43pm UTC](https://fenicsproject.discourse.group/t/temperature-dependent-material-data-from-file/4975/5 "2021-01-28T14:43:28Z")

</div>

Hmmm … since the material data change after each time step, I don’t think that storing the data once is sufficient. I just learn to understand how Function and UserExpression operates. I now switched over to:

```
...

lambdaU = Function(V)
rhoU = Function(V)
cpiU = Function(V)

lambdaU.vector()[:] = updateMatdata(lambdaX,U)
rhoU.vector()[:] = updateMatdata(rhoX,U)
cpiU.vector()[:] = updateMatdata(cpiX,U)

...

```

which avoids calling a function at every time step and I think this is at the moment ok for me. `updateMatdata` ist more or less comparable to `class userMaterial(UserExpression)`. I hope I use the Function statement the right way.