20

You should definitely check out Julia. Julia is a programming language which is similar to Python or MATLAB but utilizes a strong type-inference algorithm + JIT in order to optimize code. If types can be fully inferred in a function (which it usually can), then the code compiles in a fully static manner that matches C or Fortran performance. Translating ...


12

No, positive definiteness (and symmetry) are only precondition to using the Conjugate Gradient method. But there are plenty of other iterative methods such as MinRes and GMRES that can be used for indefinite and non-symmetric matrices.


11

What is it you want to achieve? If you want to develop methods/algorithms you might prefer a language that is flexible, and that you are familiar with. As you stated in your question, the Fortran code of your professor was hard to grasp for you, so you re-implemented it in MATLAB. This is very natural way of doing method development: first you implement ...


11

Start simple. Learn Python. I have been paid to write programms for over forty years and I have used all the languages mentioned in other answers (except Julia - I had never heard of it before now). Each language has its strengths, and most have their weaknesses. Like human languages, code is a way of expressing and framing ideas and when you know how to ...


9

Just as an aside, your github documentation is fantastic. This is just a guess from DG methods, which can have similar issues if numerical fluxes aren't chosen carefully (I figure FV methods are a subset of DG methods). If you're using interpolation from cell centers to define your fluxes, then this should be equivalent to using the average as a numerical ...


9

The finite element method is actually quite widely used in fluid flow problems, for example for the Stokes and Navier-Stokes equations. The delineation between the methods is more along the following line: The finite element method is very well suited for second order (in space) differential equations. That has something to do with the fact that trial and ...


8

The difference between finite volumes and finite differences is really more about the form of the equations solved. In typical FV methods, the conservative form is discretised in terms of integrals and fluxes, whereas FD methods generally approximate derivatives in the non-conservative form directly. Maintaining conservation and preserving physical ...


7

The basic difference is simply the meaning to be attached to the results. FDM predicts point values of any aspect of the solution. Interpolation between these values is often left to the imagination of the user. FVM predicts averages of conserved variables within specific control volumes. Hence it predicts the integrated conserved variables and can be shown ...


7

We have the following problem: $$\frac{\partial u}{\partial t}+v\color{red}{\frac{\partial u}{\partial x}}-\nu\color{blue}{\frac{\partial^2u}{\partial x^2}}=0 \tag{*}$$ The function $u$ may represent for example the concentration that propagates at velocity $v>0$ and disperses in a medium with viscosity $\nu>0$. Since only we are discussing how terms ...


7

I highly recommend to anyone, regardless of background, learning both low-level, "fast" language (C, Rust, C++, Go) high-level scripting language (Python, MATLAB, Mathematica, R, bash) As a general statement, lower-level languages far more precisely control hardware and are very efficient, while higher-level languages are easier to write/describe logic in. ...


7

Yes. The constants that appear in the interpolation estimates upon which finite element error estimates are based contain minimum and maximum angles of triangles/tetrahedra (or similar geometric measures for quadrilaterals/hexahedra). These constants are smallest whenever you have equilateral triangles or square quadrilaterals.


6

You can't. You need to approximate the integral via quadrature and or other approximations. See the following section in the document you cite.


6

The complexity of simplicial tessellations in $\mathbb{R}^{d}$, for $d>2$ is, unfortunately, not linear. In general, a simplicial tessellation $\mathcal{T}(X)$ for $n$ vertices $X\subset\mathbb{R}^{d}$ can have $\Theta(n^{\lceil d/2\rceil})$ d-simplexes in the worst case. Similar results exist for the lower dimensional faces of $\mathcal{T}$. Based on ...


6

Of course you can discretize the two equations of your system with two different methods. The challenge will simply be when the solution of one equation enters that of another. At that point, you will have to decide what the finite volume solution should be at a finite difference point should be, or how to integrate the finite difference solution (which is ...


6

The restrictive CFL of DG schemes typically comes from the combination of high order accuracy and a compact stencil (see this reference for example). The CFL depends on bounding the variational form in terms of the $L^2$ norm of the solution, which depends on derivative and traces of polynomials. Bounds for each of these quantities (using Bernstein or ...


6

The issue with finite elements in the current context is not that you have a first order differential equation, but with the kind of first order equation you have. In general, finite element methods can be made to perform quite well for first order problems such as advection equations. But that's beside the point here. To illustrate the issue, think of the ...


6

Of course, the Fourier transform is a linear operator. So, you have the kinetic energy defined as: $E(\mathbf{r}) = \frac{1}{2} \mathbf{u}(\mathbf{r}) \cdot \mathbf{u}(\mathbf{r})$. The Fourier transform of $\mathbf{u}$ and $E$ are: $$\tilde{E}(\mathbf{k}) = \int_{\Omega} E(\mathbf{r}) e^{-i \mathbf{k} \cdot \mathbf{r}} d^{3} \mathbf{r}$$ $$\tilde{\mathbf{...


6

You have some great answers already. I think there is no single answer to your question. What language(s) you choose to learn depends on what you intend to do. When I was a graduate student I too learned some Matlab, Mathematica, Maple, SAS, Stata, because my university had licenses and because in my area of research that's what others used, so I could ...


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

As Jed says, limiters are not usually an efficient approach for parabolic/elliptic problems. WENO is much more expensive than simple piecewise-polynomial interpolation, so I would first try vanilla interpolation and see if you actually have oscillations. WENO is really designed for situations in which the solution is discontinuous; yours is not. In case ...


5

The sources you are looking at are all looking at hyperbolic problems. The issues are different for elliptic problems and "limiters" are generally not the preferred tool. I outlined some of the methods and tradeoffs in this answer. As for time integration, $L$-stability is the important property to prevent bad overshoots for parabolic systems. ...


5

There are whole books written on this, but you should investigate (search the web, really) upwind diffusion and Streamline Upwind Petrov Galerkin finite element methods first. There many more methods that extend these ideas or approach the problem from a different direction.


5

The finite element method has similar problems to FD with regards to stability in solving the advection-diffusion equation, i.e. the same restrictions on Peclet number apply. One remedy, also similar to that used in FD, is to use a formulation that includes "upwinding". A nice set of lecture notes that discusses how upwinding is added to FE formulations is ...


5

If you wish to know only the change in density you are interested in the dilatation or volumetric strain of the stress tensor associated with the transformation. If you also wish to know about the anisotropy of the stretching you should look at computing the right Cauchy-Green deformation tensor, $\Delta$, and it's associated eigenvalues and eigenvectors. ...


5

The concept of $p$-nonconforming meshes can also apply to continuous FEM, not just DG. For continuous FEM, continuity is enforced by enforcing a single set of degrees of freedom on a face shared between two elements. If those elements have varying polynomial degrees, the trace space on the face must be made the same. This may be done by restricting the ...


5

"Spectral methods" usually means methods which make use of global basis functions. Fourier spectral methods use sine and cosine, and are used when you have periodic boundary conditions. Chebyshev methods use Chebyshev polynomials and are useful in non-periodic cases. These two methods are used in DNS, see e.g., hit3d which uses fourier and periodic bc, and ...


5

Taking the average will certainly work, but is recommended to take the harmonic mean to account for the changing conductivity: $k = \Big( \frac{1-f}{k_i} + \frac{f}{k_{i+1}} \Big)^{-1}$ If you choose your faces directly midway between your two cells, f = 0.5: $ k = \frac{2k_ik_{i+1}}{k_i + k_{i+1}} $ It can be derived by looking at the physical flux, ...


5

You can't. Nonlinear systems of equations are in general not solvable exactly. What you need to do is to use a method to solve nonlinear systems, of which there are of course quite a lot: A simple approach would be to use $D(U)\approx D(U^{n-1})$, where $D^{n-1}$ is the solution of the previous time step. A possibly smarter approach would be to use $D(U)\...


5

Given what you've said, I would learn C++. For one, it allows you to use MPI and lots of libraries for FEM, such as Deal.ii (which all members of this forum are obligated to mention as per our contract with Prof. Wolfgang /s). Also, if you're using C++ I imagine it would be easier to pair with openfoam (I don't know for sure as I work within my own Fortran ...


4

I have 1D finite-volume code written in python for a cell-centred mesh, First generate a sequence which is the location of the faces, $ \{ x_{j-1/2} \}$. For example, for uniform spacing over the domain [0,1] this is a simple as, a = 0 b = 1 J = 50 faces = numpy.linspace(a, b, J} After you have the faces the rest is easy because the faces uniquely ...


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