### compute expression over all dofs

90

views

0

Hi there,

I have solved a PDE and stored the solution in .h5 format file. Now i want to retrieve each time step , compute an expression over all dofs and then project the computed expression on a function space. It can be represented as following

Step 1. Solve PDE and store solution as .h5 format

Step 2. Compute expression using retrieved time step of .h5 file

Step 3. project computed expression on a function space of the same mesh as that of PDE of step1

For param = 1 the process proceeds in a satisfactory manner however for param = 70000 or even much less number like 100 i start to get very large numbers which ultimately gives nan or inf ... Can any one guide what could be the mistake ?

A MWE is could be ..

I have solved a PDE and stored the solution in .h5 format file. Now i want to retrieve each time step , compute an expression over all dofs and then project the computed expression on a function space. It can be represented as following

Step 1. Solve PDE and store solution as .h5 format

Step 2. Compute expression using retrieved time step of .h5 file

Step 3. project computed expression on a function space of the same mesh as that of PDE of step1

For param = 1 the process proceeds in a satisfactory manner however for param = 70000 or even much less number like 100 i start to get very large numbers which ultimately gives nan or inf ... Can any one guide what could be the mistake ?

A MWE is could be ..

```
from dolfin import *
mesh = Mesh('boxmesh.xml')
Q = FunctionSpace(mesh, "CG", 1)
F = Function(Q)
T = Function(Q)
T_pre = Function(Q)
T = interpolate(Expression('0',degree=1),Q)
T_pre = interpolate(Expression('0',degree=1),Q)
hdf = HDF5File(mesh.mpi_comm(), "v.h5", "r")
attr = hdf.attributes("function")
nsteps = attr['count']
for i in range(nsteps):
dataset = "function/vector_%d"%i
attr = hdf.attributes(dataset)
print 'Retrieving time step:', attr['timestamp']
hdf.read(F, dataset)
myExpr = Expression('param*((1/(1+0.05*(0.3+(1-0.1)*exp(-exp(0.05*(a - 0))))))*(b+0.05*(0.3+(1-0.1)*exp(-exp(0.05*(a - 0))))*0.005*(a-(-94.75))))',a=F,b=T_pre,param=70000,degree=1)
LagrangeInterpolator().interpolate(T, myExpr)
T_pre.assign(T)
T_out = T.compute_vertex_values(mesh)
print(T_out)
```

Community: FEniCS Project

Please login to add an answer/comment or follow this question.

just a side note.. under what circumstance would increasing degree of interpolation help me to achieve better (more accurate ) results ?

I have used degree one because the v.h5 comes from a function space which was defined over lagrange element with degree = 1. Will increasing degree of expression (in this case) lead to more accurate results ?