solving all time steps at once

Hi, the problem here is that velocity[0] etc takes out the first component
of the velocity and this is not what you
want.

Consider, in your example,
M_Mat = inner(velocity[i], phi[i])*dx
M_Mat2 = inner(velocity, phi)*dx

which should give the same, but don't:
print M_mat
print M_mat2
give
{ ((v_{-2})[0]) * ((v_{-1})[0]) } * dx0
{ ([(v_{-2})[0], (v_{-2})[1]]) : ([(v_{-1})[0], (v_{-1})[1]]) } * dx0

print out velocity also and you see what the different symbols are.

Kent

On 4 June 2012 19:00, Melanie Jahny <email address hidden>wrote:

> New question #199352 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/199352
>
> I'd like to solve the instationary stokes-equation. I solved it using a
> loop over
> each time step and it works. Now, I want to solve all time steps at once
> (there are other equations also involved later on, but this is my first
> step in
> solving my whole problem).
> The problem is that the matrix is always singular and I don't find my
> mistake.
> Is it possible to compute more than one time step at once?
> Or can anybody help me where my mistake is?
>
> Here is the code:
>
> from dolfin import *
> mesh = Rectangle(0,0,1,1,100,100)
> mesh.order()
>
> num_time_steps = 1
> dt = 0.1
>
> V = MixedFunctionSpace([VectorFunctionSpace(mesh, "CG",
> 2)]*(num_time_steps))
> W = MixedFunctionSpace([FunctionSpace(mesh, "CG", 1)]*(num_time_steps))
>
> VW = V*W
> solution = Function(VW)
> (velocity, pressure) = TrialFunctions(VW)
> (phi, psi) = TestFunctions(VW)
> vel_prev = Function(V)
> vel_prev.vector()[:] = 1.0 # just for testing
>
> # -----
> # just written here for one time step
> # later, I want to loop over more time steps here
> # -----
> i = 0
> a_st = inner(grad(phi[i]), grad(velocity[i]))*dx
> a_p = div(phi[i])*pressure[i]*dx
> B_div = psi[i]*div(velocity[i])*dx
> M_Mat = inner(velocity[i], phi[i])*dx
> M_Mat_Prev = inner(vel_prev[0], phi[i])*dx
>
> A_Mat = dt*(a_st)
> B_Mat = dt*a_p
>
> F = M_Mat - M_Mat_Prev + A_Mat - B_Mat + B_div
> a_st = lhs(F)
> L_st = rhs(F)
>
> # ---
> # boundary conditions
> # ---
> def boundary_noslip(x, on_boundary):
> return x[1]>=1.0-DOLFIN_EPS and on_boundary
>
> noslip_const = Constant((0.0, 0.0))
> noslip_top = DirichletBC(VW.sub(0), noslip_const, boundary_noslip)
>
>
> def boundary_noslip_down(x, on_boundary):
> return (x[1]<= 0.0+DOLFIN_EPS and on_boundary)
>
> noslip_down = DirichletBC(VW.sub(0), noslip_const, boundary_noslip_down)
>
> # inflow from left side
> u_in = Constant((1.0, 0.0))
>
> def boundary_noslip_inL(x, on_boundary):
> return x[0] <= 0.0+DOLFIN_EPS and on_boundary
>
> inflow_left = DirichletBC(VW.sub(0), u_in, boundary_noslip_inL)
> bcs = [noslip_top, noslip_down, inflow_left]
>
> A_st = assemble(a_st)
> b_st = assemble(L_st)
> for condition in bcs: condition.apply(A_st, b_st)
>
> solve(A_st, solution.vector(), b_st)
>
> This is the error:
> UMFPACK problem related to call to numeric
> *** Warning: UMFPACK reports that the matrix being solved is singular.
> UMFPACK problem related to call to solve
> *** Warning: UMFPACK reports that the matrix being solved is singular.
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Thank you Kent for this good hint. I see why it doesn't work.
But I don't know how to access the single components of
TrialFunctions and TestFunctions.
I tried with the "split"- function, but it doesn't work for
them.
What do I have to write in the bilinear forms in my case?
Thanks for your help.

I think I need this formulation of the mixed space:

V = VectorFunctionSpace(mesh, "CG", 2)
W = FunctionSpace(mesh, "CG", 1)

VW = MixedFunctionSpace([V, W])
Vtime = MixedFunctionSpace([VW]*(num_time_steps))

Vsub0 = Vtime.sub(0)

(vel0, pres0) = TrialFunctions(Vsub0)
(phi0, psi0) = TestFunctions(Vsub0)
velprev2 = Function(V)

But when I try to assemble the bilinear forms I get the following error:

a_st = inner(grad(phi0), grad(vel0))*dx
a_p = div(phi0)*pres0*dx
B_div = psi0*div(vel0)*dx
M_Mat = inner(vel0, phi0)*dx
M_Mat_Prev = inner(velprev, phi0)*dx

A_Mat = dt*(a_st)
B_Mat = dt*a_p

F = M_Mat - M_Mat_Prev + A_Mat - B_Mat + B_div
a_st = lhs(F)
L_st = rhs(F)

A_st = assemble(a_st)
b_st = assemble(L_st)

  File "problem.py", line 64, in <module>
    A_st = assemble(a_st)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.py", line 182, in assemble
    finalize_tensor)
  File "/usr/lib/python2.7/dist-packages/dolfin/cpp.py", line 19644, in assemble
    return _cpp.assemble(*args)
RuntimeError: *** Error: Cannot determine ownership range for sub-dofmaps.

I'm not quite sure what the best way to do this to make sure the numbering
of the unknowns is correct
and easily accessible. I'll check.
An alternative strategy could be to use block matrices in cbc.block (
https://launchpad.net/cbc.block)

Kent

On 5 June 2012 13:00, Melanie Jahny <email address hidden>wrote:

> Question #199352 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/199352
>
> Melanie Jahny gave more information on the question:
> I think I need this formulation of the mixed space:
>
> V = VectorFunctionSpace(mesh, "CG", 2)
> W = FunctionSpace(mesh, "CG", 1)
>
> VW = MixedFunctionSpace([V, W])
> Vtime = MixedFunctionSpace([VW]*(num_time_steps))
>
> Vsub0 = Vtime.sub(0)
>
> (vel0, pres0) = TrialFunctions(Vsub0)
> (phi0, psi0) = TestFunctions(Vsub0)
> velprev2 = Function(V)
>
> But when I try to assemble the bilinear forms I get the following error:
>
> a_st = inner(grad(phi0), grad(vel0))*dx
> a_p = div(phi0)*pres0*dx
> B_div = psi0*div(vel0)*dx
> M_Mat = inner(vel0, phi0)*dx
> M_Mat_Prev = inner(velprev, phi0)*dx
>
> A_Mat = dt*(a_st)
> B_Mat = dt*a_p
>
> F = M_Mat - M_Mat_Prev + A_Mat - B_Mat + B_div
> a_st = lhs(F)
> L_st = rhs(F)
>
> A_st = assemble(a_st)
> b_st = assemble(L_st)
>
> File "problem.py", line 64, in <module>
> A_st = assemble(a_st)
> File "/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.py", line
> 182, in assemble
> finalize_tensor)
> File "/usr/lib/python2.7/dist-packages/dolfin/cpp.py", line 19644, in
> assemble
> return _cpp.assemble(*args)
> RuntimeError: *** Error: Cannot determine ownership range for sub-dofmaps.
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

I tried the following code and it worked:

num_time_steps = 2
V = VectorFunctionSpace(mesh, "CG", 2)
W = FunctionSpace(mesh, "CG", 1)
VW = MixedFunctionSpace([V, W])
Vtime = MixedFunctionSpace([VW]*(num_time_steps))

(sol1, sol2) = TrialFunctions(Vtime)
(test1, test2) = TestFunctions(Vtime)

velx = sol1[0]
vely = sol1[1]
pres = sol1[2]
phix = test1[0]
phiy = test1[1]
psi = test1[2]

a_stx = inner(grad(phix), grad(velx))*dx
a_sty = inner(grad(phiy), grad(vely))*dx
a_st = a_stx + a_sty

and so on...

I still have no solution how I can do the step
(sol1, sol2) = TrialFunctions(Vtime)
(test1, test2) = TestFunctions(Vtime)
automatically for a variable number of time steps,
but so far I'm happy that it works.

I succeeded with indexing the objects
It works,
thanks for your help.