The code doesn’t give an error and it doesn’t evolve with time.
Below there is both a
PointSource and a velocity vector. Before I add PointSource, The code was running good. I am ran out of ideas.
There might be a problem due to solver, or some other thing that I couldn’t realize.
Also, where can I learn how to pick an appropriate solver for an equation? I look for “Numerical Linear Algebra with Applications” sources (books lecture notes etc.) but if you would suggest one that would be great?
What I did was the following:
#BUILD2 - pointsource # PART - 1/2 # assemble before applying pointsource A, b = assemble_system(a, L, bcs) delta = PointSource(V, Point(4., 4.), -1) delta.apply(b) # PART - 2/2 # if you define solution as uh = Function(V) # Solve PDE system at time step solve(A, uh.vector(), b)
before I was using
# Solve variational problem for time step solve(a == L, u, bcs)
The whole code:
# Import Packages from dolfin import * from mshr import * import time start_time = time.time() # Time Stepping tmax = 60 num_steps = 200 # number of time steps dt = tmax/num_steps # time step size # meshs mesh = RectangleMesh(Point(0.0, 0.0), Point(8.0, 8.0), 160, 160, "right") #1D -> 2D P1 = FiniteElement('CG', triangle, 1) # Function space & trial function - concentration V = FunctionSpace(mesh, P1) u = TrialFunction(V) v = TestFunction(V) # BUILD1- beta p1 # Define function space for velocity W = VectorFunctionSpace(mesh, 'P', 2) w = Function(W) # velocity beta = project(Expression(('-2','0'), degree=1), W) # Initial condition #1D -> 2D u_0 = Expression(('3.0'), degree=1) u_n = interpolate(u_0, V) # Boundary #1D -> 2D def boundary(x): return near(x, 0.0, DOLFIN_EPS) or near(x, 8.0, DOLFIN_EPS) or near(x, 0.0, DOLFIN_EPS) or near(x, 8.0, DOLFIN_EPS) # Condition u10 = Expression('1.0', degree=1) bcs = DirichletBC(V, u10, boundary) # PDE - Crank Nicolson F = (u - u_n)/dt*v*dx + 0.5*(Constant(0.01)*dot(grad(u), grad(v))*dx + dot(beta, grad(u))*v*dx + Constant(2.0)*v*dx + Constant(0.01)*dot(grad(u_n), grad(v))*dx + dot(beta, grad(u_n))*v*dx + Constant(2.0)*v*dx) a = lhs(F) L = rhs(F) # assemble before applying pointsource A, b = assemble_system(a, L, bcs) delta = PointSource(V, Point(4., 4.), +4) delta.apply(b) # Solution u = Function(V) # Create VTK files for visualization output vtkfile = File('bbb/bbb.pvd') # Time-stepping t = 0 for n in range(num_steps): # Update current time t += dt # Solve variational problem for time step solve(A, u.vector(), b, 'bicgstab', 'petsc_amg') # Save solution to file (VTK) vtkfile << (u, t) # Update previous solution u_n.assign(u) # Print update print ("\nTime " + str(t) + " of " + str(tmax) + "\n")
Any suggestion would be appreciated,