Using output from other simulation

I/O
1

Hello to everybody,
I want to know if it is possible to use an output of a simulation inside another simulation.
My problem is this: I’m trying to compute the symmetric part of the gradient of the velocity field u, coming from solving the stokes equation inside an obstacle with certain boundary conditions. Then I want to solve the same equation( then use the same code) but with different boundary conditions. This lead to obtain a new field, calling u_2, different from the previously u. Finally i want to compute the symgrad of u_2 and computing the volume integral between this two symgrad of u and u_2.
Is there any way to do this?

Many thanks

2

See: Loading xdmf data back in - #4 by dokken

3

Thanks for the help. I’m trying to see that question but it’s different from me. What I want is to save a function (f_32_0) called u ( which is the velocity field given by stokes system ) inside the code. But I want to save exactly with the same format (f_32_0), not in xdmf or txt or something else, because doing this i’m not able to compute the symmetric grad of f_32_0.

4

Could you provide a minimal example, as I do not understand what you are trying to do, or where you are obtaining an error.

5 
from dolfin import *
import numpy as np

mesh = UnitCubeMesh(10, 10, 10)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
B = FiniteElement("Bubble",   mesh.ufl_cell(), mesh.topology().dim() + 1)
V = VectorElement(NodalEnrichedElement(P1, B))
Q = P1
W = FunctionSpace(mesh, V * Q)


class Everywhere(SubDomain):
    def inside(self, x, on_boundary):
        return True


w = Function(W)

x = SpatialCoordinate(mesh)
W0 = W.sub(0).collapse()
u_ = project(as_vector((x[0], 2*x[1], 3*x[2]**2)), W0)
fa = FunctionAssigner(W.sub(0), W0)
fa.assign(w.sub(0), u_)

u, p = w.split()


def symgrad(u_1):
    return sym(grad(u_1))
def symgrad(u_2):
    return sym(grad(u_2))

K[1,2] = assemble(2*mu*inner(symgrad(u_1), symgrad(u_2))*dx)
print("K_12=", K[1,2])

The first time I run the code, I set a velocity along z-axis at the sphere boundaries (which is my internal boundary) and I get solution (u_1,p_1). Then I want to impose the velocity along y-axis (this means run another simulation) and I get different results calling (u_2,p_2) . Now I have two different solution and I want to compute the symgrad of both of them, so I’m looking for a method which save my previously solution doesn’t changing its structure (function). Whit a MWE is not clear what is u for me, so I specify it is the velocity field in R^3 around the obstacle.

6

Please make your code a minimal example, using built in meshes. Remove all unnecessary code to reproduce the error you are getting.

As you have not supplied how you are currently trying to use the output w (or u, p) as input to your problem, it is not clear what error you are getting. As you have defined all you boundary conditions with constants, you cannot reuse the same code to apply bcs with a function.

1 Like
7

I’ve updated the question with MWE

8

Consider the minimal working code, where you do your operation for two different inputs, and compute the combined assembly:

from dolfin import *
import numpy as np

class Everywhere(SubDomain):
    def inside(self, x, on_boundary):
        return True


def simulation(W: FunctionSpace, boundary_condition):
  mesh = W.mesh
  w = Function(W)

  W0 = W.sub(0).collapse()
  u_ = project(boundary_condition, W0)
  fa = FunctionAssigner(W.sub(0), W0)
  fa.assign(w.sub(0), u_)

  u, p = w.split()
  return u, p


def symgrad(u):
      return sym(grad(u))
 

mesh = UnitCubeMesh(10, 10, 10)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
B = FiniteElement("Bubble",   mesh.ufl_cell(), mesh.topology().dim() + 1)
V = VectorElement(NodalEnrichedElement(P1, B))
Q = P1
W = FunctionSpace(mesh, V * Q)
x = SpatialCoordinate(mesh)

u1, p1 = simulation(W, as_vector((x[0], 2*x[1], 3*x[2]**2)))

u2, p2 = simulation(W, as_vector((x[0], 2, 0)))

mu = 1
K = assemble(2*mu*inner(symgrad(u1), symgrad(u2))*dx)
print("K_12=", K)
9

sorry but I’m not understanding.

you write

u1, p1 = simulation(W, as_vector((x[0], 2*x[1], 3*x[2]**2)))

but I don’t really know why you should choose (x[0], 2x[1], 3x[2]**2). My u_1 must be equal to the u field solution obtained in the previously simulation, not be the same as any generic vector

10

The example was to illustrate how you should structure your code. You can easily replace it with:

u1, p1 = simulation(W, as_vector((x[0], 2*x[1], 3*x[2]**2)))

u2, p2 = simulation(W, u1)
11

I tried this approach :

# Write `u` to a file:
u1 = HDF5File(MPI.comm_world,"u_1.h5","w")
u1.write(u,"u_1")
u1.close()

# Read the contents of the file back into a new function, `f1`:
f1 = Function(W.sub(0).collapse())
u1 = HDF5File(MPI.comm_world,"u_1.h5","r")
u1.read(f1, "u_1")
u1.close()

# Save solution in paraview format
f1file_pvd = File("Sphere_velocity_1.pvd")
f1file_pvd << f1

Hence, what i’m expecting is to find a dolfin object u_1.h5 in the python project folder which correspond to the solution of stokes. Then I can change my boundary conditions and obtain a new field u_2, through f1 I can take the previous solution ( with different boundary conditions ) and do computations on both fields (f1, u_2)

Is this right? thanks for the patience.

12

You need to sit down and experiment, and see if you can either:
a) produce a minimal example that shows what you want to do and fails to do so.
b) find a solution that solves your problem.

2 Likes