Kyle, you will get faster answers if you post the error messages you get in addition to the program. :)
Martin On 3 May 2011 18:27, Kyle <question155...@answers.launchpad.net> wrote: > Question #155708 on DOLFIN changed: > https://answers.launchpad.net/dolfin/+question/155708 > > Kyle requested more information: > After reading Johan's post I was able to get past the nonlinear problem > but now my program crashes. Could anyone help me with this? > > Here is my program: > > #!/usr/bin/python > > from dolfin import * > import numpy, sys > > # MESHING > mesh = Rectangle(0.0, -0.5, 1.0, 0.5, 10, 10, "right/left") > > # FUNCTION SPACE > scalar = FunctionSpace(mesh, "CG", 1) #Pressure > vector = VectorFunctionSpace(mesh, "CG", 2) #Velocity > system = vector * scalar #Mixed Function Space > > # BOUNDARIES > right = compile_subdomains('x[0] == 1.0') > left = compile_subdomains('x[0] == 0.0') > top = compile_subdomains('x[1] == 0.5') > bottom = compile_subdomains('x[1] == -0.5') > > # INFLOW VELOCITY BC FOR TOP > bc1 = DirichletBC(system.sub(0), Constant((1.0, 0.0)), top) > > # NO-SLIP BC FOR TOP-BOTTOM > noslip = Constant((0.0, 0.0)) > bc0 = DirichletBC(system.sub(0), noslip, left) > bc3 = DirichletBC(system.sub(0), noslip, bottom) > bc2 = DirichletBC(system.sub(0), noslip, right) > > # STORAGE OF BCs > bcs = [bc1, bc2, bc0, bc3] > > # VARIATIONAL PROBLEM > ts = TestFunction(system) > ds = TrialFunction(system) > s = Function(system) > > v, q = split(ts) > dv, dq = split(ds) > u , p = split(s) > > f = Constant((0, 0)) > nu = 0.001 > > #F = (nu*inner(grad(v), grad(u)) - div(v)*p + q*div(u) - inner(v,f))*dx # > Stokes > F = (nu*inner(grad(v), grad(u)) - div(v)*p + q*div(u) + > inner(dot(u,grad(u)),v) - inner(v,f))*dx # Navier-Stokes > J = derivative(F,s,ds) > > # SET UP PDE > problem = VariationalProblem(F, J, bcs) > #problem.is_nonlinear=True > > # SOLVE PDE > (U, P) = problem.solve().split() > > # FILES FOR PARAVIEW > ufile_pvd = File("velocity.pvd") > ufile_pvd << U > pfile_pvd = File("pressure.pvd") > pfile_pvd << P > > # PLOTS > plot(U) > plot(P) > interactive() > > -- > You received this question notification because you are a member of > DOLFIN Team, which is an answer contact for DOLFIN. > > _______________________________________________ > Mailing list: https://launchpad.net/~dolfin > Post to : dolfin@lists.launchpad.net > Unsubscribe : https://launchpad.net/~dolfin > More help : https://help.launchpad.net/ListHelp > _______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : dolfin@lists.launchpad.net Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
