cbc team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

Interesting, how does it compare to the existing solver?

I expect it to converge slowly (in particular the solution of the
linear saddle-point problem in each iteration).

--
Anders

--- Begin Message ---

• To: Anders Logg <logg@xxxxxxxxx>
• From: noreply@xxxxxxxxxxxxx
• Date: Fri, 19 Mar 2010 23:22:15 -0000
• Delivery-date: Sat, 20 Mar 2010 00:22:19 +0100
• Reply-to: noreply@xxxxxxxxxxxxx
• Sender: bounces@xxxxxxxxxxxxx

------------------------------------------------------------
revno: 85
committer: Harish Narayanan <hnarayanan@xxxxxxxxx>
branch nick: cbc.solve
timestamp: Sat 2010-03-20 00:20:16 +0100
message:
  Added a Newton solver for the Navier-Stokes equations to the sandbox
added:
  sandbox/cg1/trial2.py


--
lp:cbc.solve
https://code.launchpad.net/~cbc-core/cbc.solve/main

Your team CBC Core Team is subscribed to branch lp:cbc.solve.
To unsubscribe from this branch go to https://code.launchpad.net/~cbc-core/cbc.solve/main/+edit-subscription.

=== added file 'sandbox/cg1/trial2.py'
--- sandbox/cg1/trial2.py	1970-01-01 00:00:00 +0000
+++ sandbox/cg1/trial2.py	2010-03-19 23:20:16 +0000
@@ -0,0 +1,61 @@
+from dolfin import *
+
+parameters["form_compiler"]["cpp_optimize"] = True
+
+mesh = UnitSquare(32, 32)
+n = FacetNormal(mesh)
+vector = VectorFunctionSpace(mesh, "CG", 1)
+scalar = FunctionSpace(mesh, "CG", 1)
+
+b = Expression(("0.0", "0.0"))
+pbar = Expression(("-1 * (1 - x[0])", "0.0"))
+
+mixed_element = MixedFunctionSpace([vector, scalar])
+V = TestFunction(mixed_element)
+dU = TrialFunction(mixed_element)
+U = Function(mixed_element)
+U0 = Function(mixed_element)
+
+xi, eta = split(V)
+v, p = split(U)
+v0, p0 = split(U0)
+
+noslip = Expression(("0.0", "0.0"))
+boundary = compile_subdomains("x[1] == 0 || x[1] == 1.0")
+bc = DirichletBC(mixed_element.sub(0), noslip, boundary)
+
+v_mid = 0.5*(v0 + v)
+p_mid = 0.5*(p0 + p)
+
+mu = Constant(1.0)
+sigma_mid = 2*mu*sym(grad(v_mid)) - p_mid*Identity(p_mid.cell().d)
+
+rho = Constant(1.0)
+dt = Constant(0.01)
+
+L = rho*inner(v - v0, xi)*dx + dt*rho*inner(grad(v_mid)*v_mid, xi)*dx \
+    + dt*inner(sigma_mid, sym(grad(xi)))*dx - dt*inner(b, xi)*dx \
+    + dt*inner(pbar, xi)*ds - dt*inner(mu*grad(v_mid).T*n, xi)*ds \
+    + dt*div(v_mid)*eta*dx
+a = derivative(L, U, dU)
+
+t = 0.0
+T = 2.0
+
+problem = VariationalProblem(a, L, bc, nonlinear = True)
+problem.parameters["newton_solver"]["absolute_tolerance"] = 1e-14
+problem.parameters["newton_solver"]["relative_tolerance"] = 1e-12
+problem.parameters["newton_solver"]["maximum_iterations"] = 50
+
+
+file = File("velocity.pvd")
+
+while t < T:
+
+    t = t + float(dt)
+
+    problem.solve(U)
+    v, p = U.split()
+    file << v
+
+    U0.assign(U)

--- End Message ---

Attachment: signature.asc
Description: Digital signature