17

If you're doing celestial mechanics over long time scales, using a classical Runge-Kutta integrator will not preserve energy. In that case, using a symplectic integrator would probably be better. Boost.odeint also implements a 4th-order symplectic Runge-Kutta scheme that would work better for long time intervals. GSL does not implement any symplectic methods,...


9

FTensor is a lightweight, header only, fully templated library that includes ergonomic summation notation. It has been tested extensively in 2, 3, and 4 dimensions, but should work fine for any number of dimensions.


9

Both GNU Scientific Library (GSL) (C) and Boost Odeint (C++) feature 8th order Runge-Kutta methods. Both are opensource, and under linux and mac they should be directly available from the package manager. Under windows, it will probably be easier for you to use Boost rather than GSL. GSL is published under the GPL license, and Boost Odeint under the ...


7

Prior answers to this question have covered most of the salient points, but I want to add one comment with respect to this: does MKL have the upper hand for some tasks? The MKL team is in a unique position to know about future Intel instruction sets and their implementations in specific processors. Furthermore, they have access to proprietary processor ...


7

While there are some books on the implementation of finite element methods (e.g., Mark Gockenbach's Understanding and Implementing the Finite Element Method and Yair Shapiro's Solving PDEs in C++: Numerical Methods in a Unified Object-Oriented Approach), I think it would indeed be most useful to look at (and compare!) actual large-scale libraries: There are ...


7

Almost everything you can build and install in your own space. With GNU autotools, you can do something like ./configure --prefix=/path/to/your/work/space ... and then follow the usual compilation instructions. Things based on CMake and Scons have similar facilities.


7

You can calculate GFLOP rates this way, but the numbers are pretty meaningless on today's hardware: Floating point operations require a variable number of clock cycles. An addition is generally cheaper than a multiplication, but each generally takes more than one clock cycle of the 2.8 billion cycles you quite. When you have hyperthreading, you have two ...


6

You can configure with separate libraries (--with-single-library=0) and link only the ones you need (e.g., -lpetscmat -lpetscvec -lpetscsys), but this is generally a waste of effort. If you use static libraries, then only the parts you reference go into the binary (if you're trying to squeeze the last megabyte out of a memory-constrained environment). PETSc ...


6

For what its worth, Eigen does have a Tensor class as an unsupported module. http://eigen.tuxfamily.org/dox-devel/unsupported/group_CXX11_Tensor__Module.html I haven't used it myself so can't say more about it. The Armadillo class library has a 3rd-order tensor class. http://arma.sourceforge.net/ I haven't used the tensor capabilities of Armadillo ...


6

I think this new taco lib is really good too. The Tensor Algebra Compiler (taco) is a C++ library that computes tensor algebra expressions on sparse and dense tensors. It uses novel compiler techniques to get performance competitive with hand-optimized kernels in widely used libraries for both sparse tensor algebra and sparse linear algebra. You can use ...


6

XTensor is a modern approach and is getting more and more popular. https://github.com/QuantStack/xtensor


6

The SINTEF Matlab Reservoir Simulation Toolbox includes a GPL-licensed AD library. The usage is mostly geared towards numerical applications in subsurface flow, but the library itself is usable for more general purposes. Here is a basic runthrough of your example as you would run it from the base directory of MRST: startup; % Load ad based module ...


6

Of course it makes sense to use the GSL (or another library for that matter) for several reasons: Don't reinvent the wheel. The work has been done, you can spend your time on more useful things. If you do decide to implement these basic things yourself, the risk that your code will probably contain some bugs and will be slower, less memory efficient etc ...


6

I would like to hear comments from users that have some practical models (e.g. black-box hyperparameter optimization) which are still needed to be solved acceptably - whether this method works or not for their models, possibly with the description of the model. Looks like you want somebody to invest what may be considerable time and energy in trying out ...


6

Optim.jl from Julia will work with the number types that you give it, so if you make it use BigFloats then it'll do that. Local derivative based, derivative-free, global, and integrates with automatic differentiation. From Julia, it's just: using Optim rosenbrock(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2 result = optimize(rosenbrock, big.(zeros(2)), ...


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

Two months after posting my question I realized that I still want to try working with deal.II, and implemented the tet->hex converter to work with complex geometry. So if someone is interested, welcome to tethex wiki page.


5

I recommend looking at different library docs: http://fenicsproject.org/book/ http://www.dealii.org/developer/


5

Ask your remote host to install what you need. We do this all the time for folks where I work. Typically they can help you out. Also, it's OK to do some of your development on your machine and then port your results to the remote machine.


5

I am a PETSc developer so take my suggestion with a grain of salt, but I would use PETSc because the problem sizes are large enough that execution overhead should be minimal, you can trivially switch between various sparse and dense solvers (1/2 sparse should be treated as dense, but it might pay off to use a sparse solver for 1/9 at your sizes), a suite of ...


5

What do people suggest for the linear system of this type? I know about trilinos, petsc, and sundials, but don't know the other alternatives or have exposure to them. Generally speaking PETSc, Trilinos, and KINSOL (from SUNDIALS) are the best-of-breed when it comes to scalability. From an ecosystem standpoint, PETSc seems most flexible, since it does not ...


5

summarizing some points: If it's a long-term integration of a non-dissapative model, a symplectic integrator is what you're looking for. Otherwise, since it's an equation of motion, Runge-Kutta Nystrom methods will be more efficient than a transformation to a first order system. There are high order RKN methods due to DP. There are some implementations, ...


5

Daniel Shapero's answer is excellent, but I felt I should add the following: Properly preconditioned iterative methods will almost certainly win here, unless your systems have a very, very special structure. There's a considerable recent literature, unfortunately mostly of the "asymptotic a priori bounds on running time" variety, on solvers for graph ...


5

You may not be able to express the weak form in deal.II as a mathematical formula, but you come pretty close. For elasticity, the bilinear form reads $$ a({\mathbf \varphi}_i, {\mathbf \varphi}_j) = \left( \lambda \nabla\cdot {\mathbf \varphi}_i, \nabla\cdot {\mathbf \varphi}_j \right)_\Omega + \left( \mu \nabla\mathbf \varphi_i, \nabla\...


5

According to [1]: "However, the IEEE-754 standard specifies nothing for elementary functions" and "Indeed, the mathematical libraries (libm) provided by operating systems do not guarantee correct rounding.". As such, according to [1], the answer is that the IEEE-754 standard does not require the exp2 function to be correctly rounded, it only recommends it. [...


4

You will likely have trouble building any package that uses the Unix/Linux build system on Windows; instead, look for packages that are built using CMake (Cross-platform Make), which is much better about being able to install packages on either platform. SparseLib++ indeed looks like it was only intended for building on Unix. To that end, I'd recommend ...


4

No. Like many other libraries, you will spend far more time trying to isolate a piece than to simply use the whole thing. Why not just install PETSc somewhere and use it as-is?


4

There is a tutorial program for this in deal.II: step-24, at http://www.dealii.org/developer/doxygen/deal.II/step_24.html (Disclaimer: I'm one of the authors of this library.)


4

I'm going to repeat Nathan's comment as an answer. You probably want LibMesh. They do have support for some Nedelec elements, though my first search found the 2D ones. If you look around or make a request, I'm sure you'll find the 3D ones, too. It has a huge element library.


4

The deal.II library (http://www.dealii.org), while written for much larger purposes, also has a sub-library of tensor classes that likely does a lot of what you want to do. In particular, it uses templates for the dimension. (Disclaimer: I am one of the principal authors of this library.)


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