AdaptiveNonlinearVariationalSolver in time stepping loop: Keeping the refined mesh

Hi Mike,

I don't have an answer for you unfortunately, just a comment: the auto-adaptive functionality in a time-loop is completely untested (at least as far as I know). In other words, you are on uncharted territory...

If you could provide a small, running example of what you are doing, I could try to comment further and perhaps provide some assistance.

Hi Marie,

Thanks for the quick reply. Here is a simple demo problem showing how I redefine the problem if the mesh has been refined in the previous time step. If I don't redefine it, then the program crashes with the message: "Dimension of function space (298) too large for application of boundary conditions to linear system (18 rows)." when it is solving the dual problem.

It is solving
du/dt = div grad u
u(0) = 1, u(1) = 0, u(x,t=0) = 0

It works fine but the code is so repetitive that I suspect there must be a better way to do it!

Thanks,
Mike

~~~~~~~~~~~~~~~
from dolfin import *

mesh = UnitSquare(5,2)
V = FunctionSpace(mesh, "Lagrange", 1) # Order 1 function space

du = TrialFunction(V)
test_u = TestFunction(V)
u = Function(V)
u0 = Function(V)

u.interpolate(Constant(0.0))

dt = 1.0
t = 0.0
tol = 1e-3

left, right, bottom, top = compile_subdomains([
 "(std::abs( x[0] ) < DOLFIN_EPS) && on_boundary",
 "(std::abs( x[0]-1) < DOLFIN_EPS) && on_boundary",
 "(std::abs( x[1] ) < DOLFIN_EPS) && on_boundary",
 "(std::abs( x[1] -1 ) < width-DOLFIN_EPS) && on_boundary"])

bcs = [DirichletBC(V, 1, left),
    DirichletBC(V, 0, right)]

F = (-(u-u0)*test_u-dt*.01*inner(grad(u),grad(test_u)))*dx
Jac = derivative(F,u,du)

problem = NonlinearVariationalProblem(F,u, bcs, Jac)
solver = NonlinearVariationalSolver(problem)
Asolver = AdaptiveNonlinearVariationalSolver(problem)
M = dot(u,u)*dx

file = File("test.pvd")

while t < (dt*20):
 print str(t)
 u0.vector()[:] = u.vector()
 #solver.solve()
 Asolver.solve(tol,M)
 t+=dt
 file << (u.leaf_node(),t)

 if u.depth>1: #If mesh was refined in previous time step
  print 'updating mesh'
  mesh = mesh.leaf_node() # save most refined mesh
  u1 = u.leaf_node().vector().array() # store refined current solution
  u01 = u.leaf_node().vector().array() # store refined past solution
  # Redefine problem based on new mesh
  V = FunctionSpace(mesh, "Lagrange", 1)
  du = TrialFunction(V)
  test_u = TestFunction(V)
  u = Function(V)
  u0 = Function(V)
  bcs = [DirichletBC(V, 1, left),
       DirichletBC(V, 0, right)]
  F = (-(u-u0)*test_u-dt*.01*inner(grad(u),grad(test_u)))*dx
  Jac = derivative(F,u,du)

  problem = NonlinearVariationalProblem(F,u, bcs, Jac)
  solver = NonlinearVariationalSolver(problem)
  Asolver = AdaptiveNonlinearVariationalSolver(problem)
  M = dot(u,u)*dx

  # Reinsert calculated solutions
  u.vector()[:] = u1
  u0.vector()[:] = u01
~~~~~~~~~~~~~~

Hi Mike,

Sorry for the lag in response -- fixing this required some surgery. If you try the development version of dolfin,
the functionality should be more as expected and hopefully give you what you are looking for too.

Note that the change has required a minor interface change: instead of calling solve with a tolerance and the goal,
the goal must now be passed to the Adaptive*Solver constructor, and solve is called with a tolerance only.

Please see below for a revised version of your original code in which no redefinition of the problem is required.

--

from dolfin import *

mesh = UnitSquare(5,2)
V = FunctionSpace(mesh, "Lagrange", 1) # Order 1 function space
du = TrialFunction(V)
test_u = TestFunction(V)
u = Function(V)
u0 = Function(V)

u.interpolate(Constant(0.0))

dt = 1.0
t = 0.0
tol = 1e-3

left, right, bottom, top = compile_subdomains([
        "(std::abs( x[0] ) < DOLFIN_EPS) && on_boundary",
        "(std::abs( x[0]-1) < DOLFIN_EPS) && on_boundary",
        "(std::abs( x[1] ) < DOLFIN_EPS) && on_boundary",
        "(std::abs( x[1] -1 ) < width-DOLFIN_EPS) && on_boundary"])

bcs = [DirichletBC(V, 1, left),
       DirichletBC(V, 0, right)]

F = (-(u-u0)*test_u-dt*.01*inner(grad(u),grad(test_u)))*dx
Jac = derivative(F,u,du)
M = dot(u,u)*dx

problem = NonlinearVariationalProblem(F, u, bcs, Jac)
solver = NonlinearVariationalSolver(problem)
Asolver = AdaptiveNonlinearVariationalSolver(problem, M)

file = File("test.pvd")

while t < (dt*20):

    # Update leaf node version of u0 (used in adaptive solve)
    u0.leaf_node().vector()[:] = u.leaf_node().vector()

    # Solve
    Asolver.solve(tol)

    # Print some stats
    info_blue("t = %g, M(u_leaf) = %g"
              % (t, assemble(dot(u.leaf_node(), u.leaf_node())*dx)))

    t+=dt
    file << (u.leaf_node(),t)

# Print summary
Asolver.summary()

Hi,
i had the same problem as Mike, but even the revised version does not work. When i want to run the above code i get the error:

    Asolver = AdaptiveNonlinearVariationalSolver(problem, M)
TypeError: __init__() takes exactly 2 arguments (3 given)

Although the documentation
http://fenicsproject.org/documentation/dolfin/dev/python/programmers-reference/fem/adaptivesolving/AdaptiveNonlinearVariationalSolver.html#dolfin.fem.adaptivesolving.AdaptiveNonlinearVariationalSolver
states, that one can give problem and the goal function...

Do i need a newer Dolfin/FEniCS version? I have installed this: http://fenicsproject.org/download/windows_details.html#windows-details

Thanks for your answer!

@Klemens: The latest binary package for Windows is version 1.0.0 and you should read the documentation for this version of DOLFIN:

http://fenicsproject.org/documentation/dolfin/1.0.0/python/programmers-reference/fem/adaptivesolving/AdaptiveNonlinearVariationalSolver.html#dolfin.fem.adaptivesolving.AdaptiveNonlinearVariationalSolver