Hot answers tagged

8

Julia has a whole ecosystem for generating sparsity patterns and doing sparse automatic differentiation in a way that mixes with scientific computing and machine learning (or scientific machine learning). Tools like SparseDiffTools.jl, ModelingToolkit.jl, and SparsityDetection.jl will do things like: Automatically find sparsity patterns from code Generate ...


5

From a performance view, you are always interested in preserving as much 'structure' in your grid as possible. Computations on a simplex- or a hexaedral mesh, where every cell looks like the next will be more performant, as you do not have to transform from local to global coordinates differently for each cell. Also, you do not have to save the cell ...


5

The short answer is that you need $$\phi_{-1} = \phi_0$$ $$\phi_N = \phi_{N-1}$$ to impose $\nabla\phi=0$. A quick check by making the following change if idx == -1: idx = 0 elif idx == N: idx = N-1 in the code, you have posted shows that the average $\phi$ remains constant up to 14 decimal places. To see why this is the correct boundary condition ...


5

As Maxim's comment points out, you ought to be able to create any solution you like, crank it through the original, continuous PDE, generate a forcing function, boundary conditions (time-dependent), and initial condition, plug those into your program, run it, and compare the answer you get to the function you started with. This is known as the Method of ...


4

The explanation in the book does not use von Neumann analysis at all but the absolute stability regions and the eigenvalues of the discrete Laplacian operator. For the result you specifically mentioned we use the fact that the maximum eigenvalues is $$ \lambda_m \approx -\frac{4}{h^2} $$ from the expression given. We then want this eigenvalue to lie inside ...


3

The choice of $k$ is restricted also by the discretization of the source term. To see it, rewrite your scheme to \begin{equation} u_m^{n+1} = \left(1 - \frac{k(1-x_m)}{h} - k(1-x_m)\right) u_m^n + \frac{k (1-x_m)}{h} u_{m+1}^n \,. \end{equation} You need $$ 1 - \frac{k(1-x_m)}{h} - k(1-x_m) \ge 0 $$ for all $x_m$. Taking $x_m=0$ (the worst case scenario) you ...


3

Most of the widely used finite element libraries are written in C++. If all you really care for -- and if all you will ever care for -- is solving an elliptic PDE on a rectangle, then it's probably not a large amount of work (a few 100 lines) to just write the finite element part yourself and use PETSc for the linear algebra. But, if you think you might ever ...


3

The notation $$\frac{\partial U}{\partial \eta}$$ means usually $$\eta \cdot \nabla U$$. This is correct even if the domain is the interval $[a,b]$. The normal vector on the interval $[a,b]$ @a is $\eta=-1$ and @b $\eta= 1$ both pointing outwards of the domain. Hence in 1D $\frac{\partial U}{\partial \eta}$ means $$\eta\cdot\nabla U=\eta \frac{dU}{dx}$$.


3

$c$ depending on time is not the issue. You will use an RK scheme which takes care of this. The issue is $c$ is discontinuous in $x$. I recommend SBP-SAT schemes for this. (1) Derive an energy equation at PDE level. (2) Search literature for SBP-SAT schemes which enforce interface conditions via SAT penalty terms, which are designed to mimic the energy ...


3

The key feature to a conservative method is simply that the changes due to the fluxes cancel out (i.e., the flux leaving one cell is entering another), so the total mass is constant. Using the form you wrote for a standard conservative method, if we sum $u$ on a grid with $M$ cells, we have: $$\sum_{m=1}^M u^{n+1}_m = \sum_{m=1}^M u^n_m - \frac{\Delta t}{\...


3

Your equation is not well posed: For general functions $f_1,f_2$, there is no function $\Phi$ so that the equation can be satisfied. For example, if you had $f_1=f_2=x$, then you are looking for a $\Phi(x,y)$ so that $$ \Phi_x = x $$ and $$ \Phi_y = y. $$ But the first of these equations imply that $$ \Phi = x^2+by+c $$ whereas the second implies that $...


3

Regarding performance, Python is definitely the bottleneck. I have experienced the same issue with a 2D Euler code I had developed, even with vectorised operations everywhere possible. It was actually even worse, as I was using solve_ivp time schemes which reallocated memory at every step... You can try and profile your code to see where the bottlenecks are. ...


3

Using the typical expansion functions (1-forms/edge-elements for E, and 2-forms/facet-elements for B) the formulations are basically the same after spatial discretization and you'd expect more or less the same accuracy. I do think they express slightly different opinions about time integration. The mixed E/B formulation nudges you in the direction of ...


3

The equations you have are what is called a "differential-algebraic equation" (DAE) because you have only time derivatives for one of the variables (namely, $u$) but not for the other (namely, $k$). The prototypical case of this kind of equation is the time dependent Stokes equations, which has time derivatives for the velocity but not the pressure....


3

As Maxim Umansky mentioned in his comment if you find the $k$ from your second equation, you would end up with this differential-integral equation: $$\frac{\partial u}{\partial t} = D \frac{\partial^{2} u}{\partial x^{2}} -5 \int_{0}^{x} u(x^{'},t) d x^{'}$$ The discretized form of this equation is: $$\frac{u_{x}^{t+\Delta t} - u_{x}^{t}}{\Delta t} = D \frac{...


3

I believe Von Neumann's stability analysis would give you the answer here. Consider the heat transfer equation: $$\frac{\partial \mathcal{T}}{\partial t} = \alpha \frac{\partial^{2} \mathcal{T}}{\partial x^{2}}$$ By using Forward Euler time integration and central difference in space discretization: $$\mathcal{T}^{t+\Delta t}_{x} = \mathcal{T}^{t}_{x} + \...


2

The potential does not pose any issues in practice. It's the Laplace term. The nonlinear term is also awkward if you have one. I've written a lengthy discussion of how one solves the (nonlinear) Schroedinger equation here. You will see that the statement "The Schrodinger equation with a zero potential is formally identical to the heat equation in the ...


2

The weak formulation is correct as stated. The space in which you are looking for solutions is $$ X = \{ v \in L_2 : \int (\partial x)^2 < \infty, v(0,y)=0, v(1,y)=0 \} $$ and this is also the space from which the test functions come. I will note that in the question, there are two other boundary conditions at the bottom and top of the box (i.e., at $y=...


2

Here is a brute force solution that would work no matter what is the discontinuity and nonlinearity in $c(x,t)$. Write your PDE as a system of two: $ \dot{y}=z\\ \dot{z}=c^2(x,t) y_{xx} $ Now, discretize it on a uniform spatial grid in x: $ \vec{x}= [x_0, x_1,..., x_{n-1}] \\ \vec{y}= [y_0, y_1,..., y_{n-1}] \\ \vec{z}= [z_0, z_1,..., z_{n-1}] \\ $ Now the ...


2

The best choice for a numerical grid is the one that will most accurately approximate the solution to your problem (without being too computationally expensive). But beyond that the specific features will depend heavily on the type of problem you are trying to solve. A grid might be aesthetically pleasing because it cleverly exploits some symmetry of the ...


2

As David said, absorbing boundary conditions won't be completely reflectionless. That said, we can reduce relfections quite a bit, which helps to avoid influence from the boundaries while the particle still travelling inside. Since this is a time dependent problem, one simple choice of boundary conditions will look like this. At the left boundary: $$\frac{\...


2

Let's take the example of the unsteady one-dimensional heat equation inside a solid on a domain $x\in[0,1]$: $$\partial_t u - D\partial_{xx} u = 0$$ with the initial profil $u(0,x) = u_0(x)$ at $t=0$, and Dirchlet boundary conditions enforcing that the wall at $x=0$ (respectively $x=1$) is at temperature $u_{L}$ (respectively $u_{R}$). If we use discretise ...


2

I've identified a few "problems": Your analytical solution isn't quite correct. The correct analytical solution to the "infinite domain" advection equation is supposed to be $u(t,x) = u_0(x-at)$. Because you have periodic boundaries, you need to properly account for this by making sure $x-at$ is properly wrapped. I replaced your ...


2

Alright I've taken a closer look a the code. One first advice that I could give you is to try to understand how the code works. It is very easy to understand that the for-loop is a time loop, i.e. the system is advanced forward one time step at each iteration of the loop. The time marching equation uses a theta-scheme for the time discretization (theta=0.5 --...


2

Let's consider the one-dimensional string first. Standard text-book physics considers the three usual boundary conditions here, namely Dirichlet (endpoints of the string are fixed), Neumann (endpoints are free) and Robin conditions (obtained e.g. when the endpoints are attached to a spring). Now, for real sonic propagation, those boundary conditions won't ...


2

I am by no means experienced with the wave equation, but I think the issue comes from the imposition of the periodic BCs. The periodic boundary conditions can be imposed by using ghost points: you do as if you were considering an extended system which, in Python terms, would have the state vector: u_extend=[u[-1], u[0], u[1], ..., u[M-1], u[M], u[0]] The ...


2

Assuming that you mistake a $t$ for an $x$ in the right hand side of the first equation, you have the following system of equations. \begin{align} \frac{\partial u}{\partial t} = D \frac{\partial^2u}{\partial x^2} + k\\ 0 = \frac{\partial k}{\partial x} - 5u\, . \end{align} To get a unique solution you would need 2 boundary conditions for $u$ and 1 boundary ...


2

You could refine your discretization and then compare the logarithm of the error ($\log |e|$) with the logarithm of the size ($\log h$) of your elements. Using a linear regression you could obtain an approximation of the order of convergence. Keep in mind that this order of convergence is asymptotic when $h \rightarrow 0$, so, for "large" $h$ you ...


2

Your nonlinear equation is a viscous Burgers' equation and could be linearized easily by using Cole-Hopf transformation. If I rewrite your equation as: $$u_{t} = 2 c \Big (\frac{\alpha}{2 c} u_{x} - \frac{u^{2}}{2}\Big)_{x}$$ And take: $$\mu = \frac{\alpha}{2c}$$ $$u = -2 \mu \frac{\phi_{x}}{\phi}$$ I would have: $$u_{x} = \frac{2 \mu}{\phi^{2}} (\phi_{x}^{2}...


1

I found a few issues with your implmentation: you should replace $t$ by $dt$ in your discrete equation (both in your code and in the question), otherwise it makes no sense ! the second part of your equation (with $a$) seems incomplete or wrong, and there's a mix between $F$ and $f$. Moreover, you flux evaluation is wrong: Fplus should read 0.5*Up**2 for ...


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