Hot answers tagged

10

I've been developing a lightweight finite element library in Python 2.7 harnessing the power of NumPy arrays and SciPy sparse matrices. The general idea is that given a mesh and a finite element, you have more-or-less one-to-one correspondence between the bilinear form and a (sparse) matrix. The user can then use the resulting matrix as he or she sees fit. ...


9

I'm a happy user of GoogleTest with a C++ MPI code in a CMake/CTest build environment: CMake automatically installs/links googletest from svn! adding tests is a one-liner! writing the tests is easy! (and google mock is very powerful!) CTest can pass command-line parameters to your tests, and exports data to CDash! This is how it works. A batch of unit-...


8

As one of the library's authors, I would of course love for deal.II to come out on top with this comparison. But I suspect it may not, and the answer lies in a factor you omit from your comparison: how long it actually takes to implement your code. Few people in academia with the skills to implement a FEM code from scratch spend more time solving PDEs than ...


6

I would recommend looking at a time-dependent example because PETSc can provide a lot more diagnostics if you formulate at that level. For example, you could use a Rosenbrock method, an additive Runge-Kutta IMEX method, or others. Some involve Newton iteration, but that is not required and is not always the best approach. To use those methods, you can ...


5

Some thoughts from someone who has worked a fair amount in compiled languages, and has done a tiny bit of FVM: Typically, if you have experience programming in C, you sketch out a high-level description (pseudocode) of what you would like to do. Then you look for libraries that might implement the data structures and capabilities you need for your high-...


5

You say that you want an MPI version. Then you need to study the literature, as the distributed memory variant of matrix-matrix product are not a simple parallellization of the sequential version. The Cannon algorithm is pretty cute if you're on a square processor grid. In each step you rotate the input matrix rows and columns, so that in the end each ...


5

Your problem is too small. You have to consider that to get good efficiency, each processor has to have enough work to offset the cost of communication. In other words, there is a threshold how many degrees of freedom there have to be below which efficiency deteriorates. I don't know where that threshold is for MUMPS. But, to give just one example, for ...


4

Heh, fun question. If you read the C standard carefully, you'll find wording like (C99, 6.5.6.8) about pointer arithmetic. If both the pointer operand and the result point to elements of the same array object, or one past the last element of the array object, the evaluation shall not produce an overflow; otherwise, the behavior is undefined. As I ...


4

If you are only interested in the smallest eigenvalue, the conjugate gradient method applied to the matrix $L$ gives you a good approximation after a reasonably small number of steps, and you won't have to solve any linear systems. The details are in Y. Saad's book on iterative methods, but here a short summary: From the coefficients that are computed in ...


4

Assuming your overall matrix is positive definite (it is definitely symmetric), then I would suggest looking into algebraic multigrid (AMG) methods as preconditioners. They compute hierarchies of sparsified matrices themselves. If you're already using PETSc, take a look at the hypre preconditioner. Using this may force you to actually multiply out the ...


4

There are several MPI-enabled software packages that use the CMake set of tools for testing. The ones that I can think of off the top of my head are Trilinos, VTK and ParaView. I would think that you don't want to assume that the executable needs to be launched with mpirun and/or mpiexec. CMake has support for specifying how to properly launch the executable ...


4

We simply roll our own code in deal.II -- in essence, we tell the framework to execute tests using mpirun -np .... We had previously just used a Makefile-based testing scheme (compile, link, execute test, then compare the output with one that had previously been saved) and you can find this here: https://svn.dealii.org/branches/releases/Branch-8-0/tests/mpi/...


4

You can either do it from the code by doing this: KSP my_solver ; // define the KSP method PC my_prec ; // define the preconditioner /* Define the matrices to be used ... */ KSPCreate(comm,&my_solver); KSPSetType(my_solver,KSPBCGS); // Sets BiCGStab as the krylov method KSPGetPC(my_solver,&my_prec); PCSetType(my_prec,PCSOR); // Sets ...


4

The PETSc team always recommends that their users control solver options from the command line. The whole package is built with the idea of extreme flexibility in composing solvers and preconditioners, and the only way to achieve that is to use the command line scheme. You could get rid of setting the KSP and PC types in your code and use ierr = KSPCreate(...


4

I like this formulation for implementing Dirichlet boundary conditions in cases where elimination of boundary dofs is not convenient. If you apply effectively the same procedure to $M$, the rows and columns corresponding to Dirichlet boundary dofs will be rows and columns of the identity. This will give you new eigenvalues equal to 1 (or $\alpha$ in the ...


4

Continuous weak form Though I think the weak form is more fundamental, suppose we start with the strong divergence-form representation for a (first- or) second-order quasilinear PDE: find the $m$-component solution $u \in R^m(\Omega)$ $$ -\nabla \cdot f_1(u,\nabla u) + f_0(u,\nabla u) = 0 $$ on the domain $\Omega \subset R^d$ with (for simplicity) ...


4

If you're using iterative methods with assembled matrices, just buy DDR channels. Don't pay attention to number of cores when loooking at the spec sheet. Within the same class of processors (e.g., a recent generation of Xeon), the achievable memory bandwidth will be proportional to the quoted peak bandwidth. Note that for very small problem sizes, you might ...


4

Unfortunately there's no tool for this. You can run each on a variety of input sizes to establish the computational complexity they appear to have, i.e. the $f$ in the $O(f(n))$ that characterizes each code. This can point you towards what underlying algorithms each is using and verify or not that something that should be $O(n)$ is actually implemented to ...


4

It is indeed possible to write an obstacle problem for an advection-diffusion equation as a variational inequality: If $a(u,v)$ is the bilinear form corresponding to your advection-diffusion equation, the corresponding obstacle problem (in your specific case) is finding $u\geq 0$ such that $$ a(u,v-u) \geq (f,v-u) \qquad \text{for all } v\geq 0.$$ In ...


4

I think you have some confusion. PETSc is not in the same league as Fenics, Libmesh, Moose etc. In fact, all of these (heavyweight) packages use PETSc for linear algebra. IMHO PETSc is as lightweight as you can get. It just requires C/Fortran compilers and Python (used only for configuration) and you can build the library in under 5 minutes on your laptop. ...


4

I can recommend nutils. nutils meets at least a few your "light-weight" requirements. it is pure python and easy to install since it only depends on standard Python libraries numpy, scipy, and matplotlib and, thus, it is well suited for interoperations. At least the developers claim that "Exposed objects are of native python type or allow for easy ...


4

Wow! A question that I can answer! I have been using PETSc for the past year and a half to solve the Navier-Stokes equations (with some hard-coded MPI). The best way to learn PETSc is to (1) read the manual so you generally know what's in PETSc, (2) decide roughly which part of PETSc you need, and then (3) begin experimenting with simple examples to ...


4

Such problems can be solved using an $LDL^T$ factorization (similar in memory and time cost to Cholesky). Your matrix is not very sparse so treating it as such may have limited benefit. I would recommend comparing sparse direct solvers (such as MUMPS and Umfpack) to the dense factorizations you've used. Complex-symmetric matrices are often difficult to ...


4

Use of Lagrange multipliers produces a saddle-point problem, $$ \begin{pmatrix} A & B^T \\ B & 0 \end{pmatrix} \begin{pmatrix} u \\ \lambda \end{pmatrix} = \begin{pmatrix} b \\ 0 \end{pmatrix} $$ As you've noticed, many preconditioners break down for this sort of system. One can use direct solvers that support pivoting, but if you want iterative ...


3

As of writing this answer (October 2014), Trilinos does not have this capability throughout its code base. There are two packages (Kokkos and Tpetra) that will provide this functionality to Trilinos, but neither are used widely throughout the remainder of the library. The Trilinos plan (as far as I understand it) foresees new packages (and some revamped old ...


3

ADI is not a very good parallel algorithm. You should seriously consider formulating the problem in 3D and solving with multigrid. You could get a start with src/ksp/ksp/examples/tutorials/ex45.c. If you insist on using ADI, you should seriously consider allocating the three matrices separately. It's more memory, but then you won't have to reassemble on ...


3

There are generally three approaches to this problem: 1/ If processor P has a deterministic algorithm by which it needs to know which matrix elements it needs, then processor 1 can replay this same algorithm to find out which elements processor P will need. It can then send these elements to processor P (and P alone, not everyone). 2/ Processor P can "...


3

It is just a projection $1 - ZZ^T$. If you are digging around in the implementation, look at KSP_PCApplyBAorAB() and related functions.


3

For linear problems of size ~1000, you can't beat a direct solver by trying to use an iterative solver. For problems of this size, it also doesn't really matter. If your problem was any bigger, I would recommend you watch/listen lectures 34 and 38 of the course I recorded this spring: http://www.math.tamu.edu/~bangerth/videos.html .


3

MatSetValuesBlocked contributes logically dense blocks. Given the pattern shown above, you should be calling it once per block-row. You're currently calling it as though the arguments were COO format. Also be aware that reading a matrix from a file is a workflow choice that guarantees that your software will never be scalable. Long-term, you should ...


Only top voted, non community-wiki answers of a minimum length are eligible