how to speed up a linear solver?

That is a huge system. I would not advise using LU to solve a system
with 12M dofs. Use an iterative solver.

--
Anders

On Tue, Oct 30, 2012 at 02:30:58AM -0000, Eric Li wrote:
> New question #212710 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/212710
>
> Hi all,
>
> Recently I realize that it is quite time consuming for solving a problem with high spatial resolution(say 1524*1024). To be more specific, I run the following SAMPLE code(solving for velocity field from stream-function) in parallel with 20 processors, and it got stuck at the point of solving linear system.
> (Process 0: Solving linear system of size 12494802 x 12494802 (PETSc LU solver, mumps)) ->> got stuck here
> I wonder if there is any way to speed it up. Please advise at any rate. Thank you
>
>
> #----Import modules-----------
> from dolfin import *
> import math
> sin=math.sin
> #-----------------------------
> set_log_level(DEBUG)
>
> X0 , D, nx , ny = 4*DOLFIN_PI , 3.0 , 1524 , 1024
>
> tol = DOLFIN_EPS
>
> mesh = Rectangle(0,-D,X0,D,nx,ny)
>
>
> file_psi = File('./psi/psi.pvd','compressed')
> file_velocity = File('./vel/velocity.pvd','compressed')
>
>
> Scalar=FunctionSpace(mesh,'CG',2)
> Vector=VectorFunctionSpace(mesh,'CG',2)
>
> #------ Define Boundaries ------------------------------------
> def top_boundary(x,on_boundary):
> return on_boundary and abs(x[1]-D) < tol
> def bot_boundary(x,on_boundary):
> return on_boundary and abs(x[1]+D) < tol
> def left_boundary(x,on_boundary):
> return on_boundary and abs(x[0]) < tol
> def right_boundary(x,on_boundary):
> return on_boundary and abs(x[0]-X0) < tol
>
>
>
> #---- Boundary Conditions for Velocity -----------------------------
>
> bc_t_vel = DirichletBC(Vector.sub(0), 0.0, top_boundary)
> bc_b_vel = DirichletBC(Vector.sub(0), 0.0, bot_boundary)
> bc_l_vel = DirichletBC(Vector.sub(0), 0.0, left_boundary)
> bc_r_vel = DirichletBC(Vector.sub(1), 0.0, right_boundary)
>
> bcs_vel =[bc_t_vel, bc_b_vel, bc_l_vel, bc_r_vel]
>
>
> psi = Expression('(1.0/3.0*x[1]*x[1]*x[1]-9.0*x[1])*(sin(x[0])-2*x[0])')
> psi = interpolate(psi, Scalar)
>
> file_psi << psi
>
>
> velocity = TrialFunction(Vector)
> test = TestFunction(Vector)
> a_vel = inner(velocity,test)*dx
> L_vel = (psi.dx(1)*test[0]-psi.dx(0)*test[1])*dx
> velocity = Function(Vector)
> solve(a_vel == L_vel, velocity, bcs_vel)
>
>
> file_velocity << velocity
>

Hi Anders,

Thank you for your answer, could you tell me which iterative solver is good in this case? I notice that some solvers may not work in parallel.

Thanks Kent-Andre Mardal, that solved my question.