>> In libMesh this is as simple as providing a residual function and then >> using the -ksp_matrix_free or -ksp_mf_operator (I think) options to PETSc. > > I understand this. But my concern is how would you compute the local > contributions to the mat-vec product ? Since I do not want to store the > global stiffness or mass matrices, the question then becomes, what is the > "overhead" in writing a user routine to perform local assemblies with local > mass, stiffness for each DOF to get the nonlinear residual functional ? And > how does this operation scale in parallel ?
Oh... That part is easy. In all the examples we show, we compute an element stiffness matrix and then assemble the global stffness matrix: K = Am(K_e) & f = Av(f_e) where Am(), Av() are notional matrix and vector assembly operators. Similarly, let the element residual be r_e = K_e u_e - f_e and the global residual is r = Av(r_e). So in the matrix free approach we just calculate the global residual from local contributions, and never form a global matrix. Now, of course global matrix assembly is expensive in terms of storage, but the vector assembly is not. For a vector of length N on P processors, each processor stores ~O(N/P) entries. And very few elements actually need to be communicated. In 2D the parallel vector summation is often driven by network latency instead of bandwidth because the data transfers are small. > I pose this question since the matrix-free algorithm will not involve any > kind of assembling routine at all. All it would need is to know the local h, > p, and basis type to compute local mass, stiffness matrices to compute the > local residual (a value in the global residual vector corresponding to the > unknown). Agreed. > Also, I perfectly understand your concern about creating a global > preconditioner matrix for this method to be successful but I believe that a > framework can be written well to handle the matrix-free algorithm in a > purely recursive fashion thereby eliminating any need for a preconditioner > matrix. Just my 2 cents. > Any pointers to benchmark studies that you have done in the past regarding > the speed and scalability for different kinds of problems (diffusion, > convection dominated) would also be very useful to me. I'll see what I can pull up. PETSc's interface to hypre for diffusion problems allows for algebraic multigrid directly, which is hard to beat. > I would be happy to contribute to a sample code regarding a matrix-free > scheme with LibMesh and PETSc which will give me an opportunity to test out > different techniques I have in mind. I was thinking about using the Laplace-Young free surface model to put together an example problem using this functionality. The LY model is particularly well-suited since it is 2D, scalar, steady, and nonlinear, and well-suited to adaptivity. I hope to get started on it Monday - I'll be on travel and that's always a good time for me to code without the daily distractions. ;-) I should have something by the end of the week. -Ben ------------------------------------------------------------------------- SF.Net email is sponsored by: Check out the new SourceForge.net Marketplace. It's the best place to buy or sell services for just about anything Open Source. http://ad.doubleclick.net/clk;164216239;13503038;w?http://sf.net/marketplace _______________________________________________ Libmesh-users mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/libmesh-users
