问:
Equivalent of multiple-scale analysis for a linear ODE
I came across the method of multiple-scale analysis and was intrigued, because I am trying to solve a linear ODE with multiple characteristic timescales. When I apply the method as described, say, here, I find that the zeroth-order solution is the exact solution! This is exactly what I don't want, because conventional methods have not been successful in integrating this equation numerically. I was very interested in the possibility of approximating the solution.

Derivation

I will briefly recount what I tried, which follows the linked document closely.

Suppose I have an operator $A = A_0 + \epsilon A_1$ in dimensionless units, where $\epsilon \ll 1$. The equation I want to solve is $\dot{y} = A y$.

Suppose I assign two timescales, $t_0 = t$ and $t_1 = \epsilon t$, as described in the document. Then, when I take the time derivative and gather terms of the same power of $\epsilon$, I get

$$\left(\frac{\partial y_0}{\partial t_0} - A_0 y_0\right) + \epsilon \left(\frac{\partial y_0}{\partial t_1} + \frac{\partial y_1}{\partial t_0} - A_1 y_0 - A_0 y_1 \right) + \mathcal{O}\!\left(\epsilon^2\right) = 0$$

The solution to the first (order unity term) is

$$ y_0 = e^{A t_0} x_0\!\left(t_1\right) $$

and, when we solve the second term with this solution substituted,

$$ y_1 = e^{A_0 t_0} x_1\!\left(t_0\right) + e^{A_0 t_0} t_0 \left(A_1 x_0\!\left(t_1\right) - \frac{\partial x_0}{\partial t_1}\right) $$

Now, solving for the parenthetical term to go to zero (because it's secular),

$$x_0 = e^{A_1 t_1} z_0$$

Substituting this back into $y_0$ from earlier, we find that

$$y_0 = e^{A_0 t_0 + A_1 t_1} z_0 = e^{\left(A_0 + \epsilon A_1\right) t} z_0$$

which is exactly the solution to the linear equation.

Other attempts

I have tried to implement a kind-of IMEX (?) scheme, where I solve

$$ e^{\epsilon A_1 h} y^{n+1} = e^{A_0 h} y^n $$

which is implicit for the fast timescale and explicit for the slow timescale. This also doesn't work, because it is only valid until the simulation time reaches the slow timescale. The next logical thing to do is to have a transition from this scheme to a fully-implicit scheme when the time is large enough, but, in practice, that is too computationally expensive.

Question

I am curious to know if I've just done this wrong, but also if there are methods in the same spirit as the multiple-scale analysis that would work for a linear ODE.


问:
Minimum of quadratic assignment (QAP) with convex objective
Suppose $A,B\succeq0$ and $C\in\mathbb R^{n\times n}$. I am hoping to solve an instance of the following optimization problem:
$$
\min_{\textrm{permutation matrices }P}
\mathrm{tr}(BP^\top AP+C^\top P)
.
$$
This objective function is convex in $P$, but the constraint forces $P$ to be a permutation matrix.  In essence, this is the quadratic assignment problem (QAP) with an additional assumption that the objective function is convex.

If $A=0$ or $B=0$, the problem would be the linear assignment problem.  If we relaxed the constraint and thought of $P$ as doubly-stochastic, we could optimize the resulting convex problem and get a permutation matrix anyway (the vertices of the Birkhoff polytope are permutations).

Is there any similar trick for the problem above?  I know QAP is NP-hard in general, but maybe it's easier if the objective function is convex?

ContentLicense: CC BY-SA 4.0


问:
Runge-Kutta for PID and system in separate calculations without filter
I need to calculate a closed-loop system in Python; specifically, obtain the PID response and then use the output to obtain the system response sample-by-sample with my own loop.

For this, I am using a fourth-order Runge-Kutta with an adaptitive step size to solve the state space representation. So far, the only way I managed to get a good aproximate response was by adding a low-pass filter with the PID:

In transfer function form, my PID without filter:

$$ G(s)_c = K_p + \displaystyle\frac{K_i}{s} + \displaystyle\frac{K_dN}{1+ \displaystyle\frac{N}{s}}$$

My PID with the filter:

$$ G(s)_c = \displaystyle\frac{1}{0.1s + 1}\left(K_p + \displaystyle\frac{K_i}{s} + \displaystyle\frac{K_dN}{1+ \displaystyle\frac{N}{s}}\right)$$

With an example system:

$$G(s) = \displaystyle\frac{1}{s^2 + s + 1}$$,

where $K_p = K_i = K_d = 1$ and $N = 100 $

Result with the filter:

and without the filter:

The expected response should resemble the first image, but the filter adds tiny changes to the response. If I calculate PID + system in one transfer function, I get the expected results and I don't need the filter (probably because the system act as a filter), but, like i said, the requirement is to make the calculation of the PID separately.

I can also obtain the expected response if I increase the relative tolerance and absolute tolerance of the RK4, but the amount of increase in the # of steps is unacceptable, from 890 to 4269. The system itself does not matter, could be anything.

Any ideas on how could I get an accurate result without the filter or a big increase in number of steps?

the complete python code:

import numpy as np
import control as ctrl
from matplotlib import pyplot as plt
import time

def norm(x):
    return np.linalg.norm(x) / x.size**0.5

def runge_kutta(ss, x, h, inputValue):
    k1 = h * (ss.A * x + ss.B * inputValue)
    k2 = h * (ss.A * (x + k1/2) + ss.B * inputValue)
    k3 = h * (ss.A * (x + k2/2) + ss.B * inputValue)
    k4 = h * (ss.A * (x+k3) + ss.B * inputValue)

    x = x + (1/6) * (k1 + k2*2 + k3*2 + k4)
    y = ss.C * x + ss.D * inputValue
    return y.item(), x

N = 100
kp = 1
ki = 1
kd = 1

pid = ctrl.tf2ss(
    ctrl.TransferFunction([1], [0.1, 1]) *
    ctrl.TransferFunction([N*kd + kp, N*kp + ki, N * ki], [1, N, 0]))

x_pidB = np.zeros_like(pid.B)
x_pidS = np.zeros_like(pid.B)

system = ctrl.tf2ss(ctrl.TransferFunction([1], [1, 1, 1]))
vstadosB = np.zeros_like(system.B)

min_step_decrease = 0.2
max_step_increase = 5
h_ant = 0.0001
rtol = 1e-3
atol = 1e-6
Time = 0
tbound = 30
sp = 1
output = [0]
Time_out = [0]
yb = 0
sf1 = 0.9
sc_t = [0]
start = time.time()
counter = 0

while Time < tbound:
    error = sp - yb
    while True:
        counter += 1
        if Time + h_ant >= tbound:
            h_ant = tbound - Time
            ypids, x_pidSn = runge_kutta(pid, x_pidS, h_ant, error)
        else:
            ypidb, x_pidBn = runge_kutta(pid, x_pidB, h_ant, error)
            ypids, x_pidSn = runge_kutta(pid, x_pidS, h_ant / 2, error)
            ypids, x_pidSn = runge_kutta(pid, x_pidSn, h_ant / 2, error)

            scale = atol + rtol * (np.abs(x_pidSn) + np.abs(x_pidBn)) / 2
            delta1 = np.abs(x_pidBn - x_pidSn)
            error_norm = norm(delta1 / scale)

            if error_norm == 0:
                h_est = h_ant * max_step_increase
            elif error_norm < 1:
                h_est = h_ant * min(max_step_increase,
                                    max(1, sf1 * error_norm**(-1 / (4+1))))
            else:
                h_ant = h_ant * max(min_step_decrease, sf1 * error_norm**(-1 / (4+1)))
                continue

        yb, vstadosB = runge_kutta(system, vstadosB, h_ant, ypids)
        break

    sc_t.append(ypids)
    output.append(yb)
    Time += h_ant
    Time_out.append(Time)

    h_ant = h_est
    x_pidB = x_pidSn
    x_pidS = x_pidSn

print(f'Number of calculations:{counter}')
print(f'Time:{time.time() - start:.3f}')

fig = plt.figure(212)
ax1 = fig.add_subplot(211)
ax1.plot(Time_out, output, label='System response')
ax1.legend()
ax1.grid()

ax2 = fig.add_subplot(212)
ax2.plot(Time_out, sc_t, 'r', label='Control signal')
ax2.legend()
ax2.grid()

plt.show()
```


问:
Inverses of many standard subspaces of one large matrix
i have a large rectangular invertible matrix M (about 5000x5000) and i have a loop in which i do the following for each iteration i (there are about 6000 iterations):

i am given a subspace S_i (which in my problem is always a subset of the standard basis; for example, "dimensions 3, and 300, and 3723";  these subspaces will consist of about 100 out of the original 5000 dimensions). I have no control over which subspace i will be given, and it is not random; some of the dimensions are more 'popular' than others and will tend to more frequently appear within S_i. In my problem, it is unlikely that S_i will happen to be an invariant subspace of M.

Compute the rectangular matrix N_i which is just the restriction of the matrix N_i to S_i (that is, start with M and then drop all of the rows and columns of the matrix except for those corresponding to the subspace; for example, drop all rows except for 3, 300, 3723 and drop all columns except for 3, 300, and 3723)

Take the inverse N_i^-1

For example, if the matrix M is:

```<code>1 2 3
4 0 6
7 8 9
</code>```

and the subspace for some particular iteration is dimensions {1,2} (with one-based indexing), then N_i is:

```<code>1 2
4 0
</code>```

and the desired inverse for this iteration is:

```<code>0    0.25
0.5 -0.125
</code>```

i feel that perhaps, since i am taking the inverse of parts of the same original matrix M over and over again, and also because many some of the rows and columns are more frequently chosen than others, there could be some shared computation between these 6000 iterations, particularly in step 3 when the inverse is computed. Is there a more efficient way to do this?


问:
On using Ritz Method to solve a Mindlin–Reissner plate
I am trying to replicate the method given in the this paper. I have written a Matlab program which determines the displacement field of Mindlin–Reissner plate theory using Ritz method. The limitation of the program is that the plate should lie on xy-plane.

Now I want to extend it to plates that are NOT aligned with any coordinate plane. One such example is given below:

My questions are:

How can I find the Jacobian matrix in this case?
$x = \sum_{i=1}^{4}x_i \times N_{i}(\xi, \eta)$ and $y = \sum_{i=1}^{4}y_i \times N_{i}(\xi, \eta)$ still valid in this case?

When integrating with respect to $z$, what will be the limits of integration? For a plate parallel to $xy-plane$, I have used $z_m - \frac{t}{2}$ to $z_m + \frac{t}{2}$, where $z_m$ is the distance of middle surface from $xy-plane$ and $t$ is the thickness of the plate.

My Matlab program can accessed from this link.


问:
Orientation of cones and transitive closures of a dmplex in Petsc
From the Petsc manual pages, I fail to understand what the orientations of cones and transitive closures mean.

In particular, how can I relate these numbers to the orientation of the length/area/volume of the corresponding mesh components ? How can I retrieve the notion of "relative orientation of a face within a cell", with which I am more familiar.

Thanks in advance


问:
Solution errors when refining a static grid: Continuous vs. step-wise refinement
Let's assume I am working on a 2-D domain on $R^2$, with my coordinates $x \in[-1,1]$, $y \in[-1,1]$ and I want to solve a popular CFD problem, like the shallow water system or the Euler system.

At $x=0$ there is some physical feature that I want to resolve really well, so I want to employ mesh refinement around the origin. I am restricted to squares as resolution element.

I know that my solution around the origin will be close to hydrostatic, but possibly evolving and rotating slowly. This could be e.g. a cyclonic system.

I would now be interested in the following:

Which static technique of mesh refinement, continuous (e.g. with $dx \propto 1-exp(-x^2)$) or step-wise (often also called nested meshes) would give me the smaller errors in my solution?

I imagine that this may have been done already in the atmospheric community, which uses global grids with local refinement, but so far I came up with only step-wise refinement techniques. Also this review paper on cubed sphere grids mentions that non-uniformity is bad for a dynamic situation, but that's it.

Any answer that would also give a citable (review) paper to answer this question would be greatly appreciated.


问:
Multipole expansion for magneticfield intesity (magnetization)
I'm using the Barnes Hut method to calculate the magnetic vector potential induced by an applied current. Given as:

$\begin{equation} A(r) = \frac{\mu_0}{4\pi} \int_V\frac{\bf{J(r')}}{|r-r'|}dV(r') \end{equation}$

In general this equation can be written as:

$\begin{equation} \Phi_i = \sum^N_{j=0} \frac{m_j}{r_{ij}} \end{equation}$

By Taylor expanding the general form $\frac{m_j}{r_{ij}}$ at a cell center I can calculate my multipole which are only dependet on $j$:

$\begin{equation}\begin{split}
\frac{m_j}{r_{ij}} & \approx \frac{m_j}{R} + m_j(x_c-x_j)\left(-\frac{x_i-x_c}{R^3}\right) + m_j(y_c-y_j)\left(-\frac{y_i-y_c}{R^3}\right) \\
& + m_j(z_c-z_j)\left(-\frac{z_i-z_c}{R^3}\right) + \frac{m_j(x_c-x_j)^2}{2}\left(\frac{3(x_i-x_c)^2}{R^5} - \frac{1}{R^3}\right) \\
& + \frac{m_j(y_c-y_j)^2}{2}\left(\frac{3(y_i-y_c)^2}{R^5} - \frac{1}{R^3}\right) + \frac{m_j(z_c-z_j)^2}{2}\left(\frac{3(z_i-z_c)^2}{R^5} - \frac{1}{R^3}\right) \\
& + \frac{m_j(x_c-x_j)(y_c-y_j)}{2}\frac{3(x_i-x_c)(y_i-y_c)}{R^5} \\
& + \frac{m_j(y_c-y_j)(z_c-z_j)}{2}\frac{3(y_i-y_c)(z_i-z_c)}{R^5} \\
& + \frac{m_j(z_c-z_j)(x_c-x_j)}{2}\frac{3(z_i-z_c)(x_i-x_c)}{R^5}
\end{split}\end{equation}$\

$\begin{equation} \Phi_i = \sum^N_{j=0} \frac{m_j}{r_{ij}} = \sum^N_{j=0} A_j B_i \end{equation}$

My question now is: can I do a multipole expansion for calculating the magnetic field intesity with consideration of the magnetization:

$\begin{equation} H(r) = H_i(r) -\frac{1}{4\pi} \int_V [\frac{M(r')}{|r-r'|^3} - \frac{3(r-r')}{|r-r'|^5}(r-r') \cdot M(r')]dV(r') \end{equation}$

Here $M$ is the magnetization of a ferromagnetic material which actually depends on the magnetic field intesity. The magnetic field intensity $H_i$ is given by Biot-Savarts law from all external electric currents.
If this is possible what would be a good strategy to get independent $A_j$ and $B_i$? I actually try to Taylor expanded the equation but the variables are not independet as you can see for the first derivative in x-direction (without ($H_i$)):

\begin{equation} H_{x_j|c} =
\newcommand{\Bold}[1]{\mathbf{#1}}\left(-\frac{6 \, M_{\mathit{xj}} {\left(x_{c} - x_{i}\right)}}{{\left({\left| -x_{c} + x_{i} \right|}^{2} + {\left| -y_{c} + y_{i} \right|}^{2} + {\left| -z_{c} + z_{i} \right|}^{2}\right)}^{\frac{5}{2}}} + \frac{15 \, {\left(M_{\mathit{xj}} {\left(x_{c} - x_{i}\right)} + M_{\mathit{yj}} {\left(y_{c} - y_{i}\right)} + M_{\mathit{zj}} {\left(z_{c} - z_{i}\right)}\right)} {\left(x_{c} - x_{i}\right)}^{2}}{{\left({\left| -x_{c} + x_{i} \right|}^{2} + {\left| -y_{c} + y_{i} \right|}^{2} + {\left| -z_{c} + z_{i} \right|}^{2}\right)}^{\frac{7}{2}}} - \frac{3 \, {\left(M_{\mathit{xj}} {\left(x_{c} - x_{i}\right)} + M_{\mathit{yj}} {\left(y_{c} - y_{i}\right)} + M_{\mathit{zj}} {\left(z_{c} - z_{i}\right)}\right)}}{{\left({\left| -x_{c} + x_{i} \right|}^{2} + {\left| -y_{c} + y_{i} \right|}^{2} + {\left| -z_{c} + z_{i} \right|}^{2}\right)}^{\frac{5}{2}}},\,-\frac{3 \, M_{\mathit{yj}} {\left(x_{c} - x_{i}\right)}}{{\left({\left| -x_{c} + x_{i} \right|}^{2} + {\left| -y_{c} + y_{i} \right|}^{2} + {\left| -z_{c} + z_{i} \right|}^{2}\right)}^{\frac{5}{2}}} - \frac{3 \, M_{\mathit{xj}} {\left(y_{c} - y_{i}\right)}}{{\left({\left| -x_{c} + x_{i} \right|}^{2} + {\left| -y_{c} + y_{i} \right|}^{2} + {\left| -z_{c} + z_{i} \right|}^{2}\right)}^{\frac{5}{2}}} + \frac{15 \, {\left(M_{\mathit{xj}} {\left(x_{c} - x_{i}\right)} + M_{\mathit{yj}} {\left(y_{c} - y_{i}\right)} + M_{\mathit{zj}} {\left(z_{c} - z_{i}\right)}\right)} {\left(x_{c} - x_{i}\right)} {\left(y_{c} - y_{i}\right)}}{{\left({\left| -x_{c} + x_{i} \right|}^{2} + {\left| -y_{c} + y_{i} \right|}^{2} + {\left| -z_{c} + z_{i} \right|}^{2}\right)}^{\frac{7}{2}}},\,-\frac{3 \, M_{\mathit{zj}} {\left(x_{c} - x_{i}\right)}}{{\left({\left| -x_{c} + x_{i} \right|}^{2} + {\left| -y_{c} + y_{i} \right|}^{2} + {\left| -z_{c} + z_{i} \right|}^{2}\right)}^{\frac{5}{2}}} - \frac{3 \, M_{\mathit{xj}} {\left(z_{c} - z_{i}\right)}}{{\left({\left| -x_{c} + x_{i} \right|}^{2} + {\left| -y_{c} + y_{i} \right|}^{2} + {\left| -z_{c} + z_{i} \right|}^{2}\right)}^{\frac{5}{2}}} + \frac{15 \, {\left(M_{\mathit{xj}} {\left(x_{c} - x_{i}\right)} + M_{\mathit{yj}} {\left(y_{c} - y_{i}\right)} + M_{\mathit{zj}} {\left(z_{c} - z_{i}\right)}\right)} {\left(x_{c} - x_{i}\right)} {\left(z_{c} - z_{i}\right)}}{{\left({\left| -x_{c} + x_{i} \right|}^{2} + {\left| -y_{c} + y_{i} \right|}^{2} + {\left| -z_{c} + z_{i} \right|}^{2}\right)}^{\frac{7}{2}}}\right) \end{equation}


问:
Numerical integration of SDE: choice of $dt$ and algorithm
I am working on the following Stochastic Differential Equation (SDE) in the Quantum Mechanics context:

$$dX_{t} = a X_{t} dt + b X_{t} dW$$

where $X_{t}$ is my stochastic varible, $dt$ is my integration time step, and $dW = N[0,1] \sqrt{dt}$ is the Wiener differential.

Due to the "quantumness" of this SDE, it's not possible to obtain an analytical solution, in contrast to the classical case (Geometric Brownian Motion) for which it is available.

So, to solve this equation numerically I have implemented Euler-Maruyama and Milstein methods. However, the choice of the time step $dt$ for the integration seems to be arbitrary since I have no comparison with an exact solution. Is there any general rule for the choice of $dt$?

I was thinking about study the limitations of the choice of $dt$ in the framework of ODEs Euler scheme for the equation without the noise term (for which I have an analytical solution) and use that to choose $dt$ for the SDE. Is it feasible?

Moreover I have to change the magnitude of the noise term to be from a perturbative term to the predominant one.
I read in another thread that it is not possible to adapt regular ODEs solver like 4th-order Runge-Kutta method to add the noise term $dW$
in an Euler/Milstein fashion, but in the case of perturbative noise the deterministic part is predominant over the noise one, and using a better approximation for the deterministic part can allow choosing bigger $dt$, saving computational time. Is still wrong to adapt 4th-order Runge-Kutta to the SDE? Is there any other method suitable for small noise?


问:
Predictor-Corrector vs. Deferred Difference Corrections
I want to use the Numerov method but keep higher-order terms from the Taylor expansions. In the literature, I found the term "Deferred Difference Corrections" for the procedure of first solving the differential equations via the standard Numerov method and afterward adding the higher-order terms. For me, this sounds just like a Predictor-Corrector method. Can someone tell me what the difference is?</p>


问:
FEniCS implementation of Maxwell equations for a dipole antenna
someone knows where I can find a FEniCS implementation of Maxwell equations for a dipole or other type of antenna? I mean a dipole antenna with an arbitrary geometry of every 'leg' in the dipole. </p>


问:
numerical instabilities in Fluid Dynamics, Finite Element Method
I'm looking for references to understand where the numerical instabilities come from in hydrodynamics in general, and notably when the Péclet number: $Pe>1$. I'm using the finite element method.

Would you have some references for that, please? Or if you could explain at least the case of the Péclet instability for example in the equation $c_t+u.\nabla c=D\nabla^2c$ in 1d, even in the finite difference method I'd be happy.

I looked at this course and it looks like it comes from numerical observation... Is there any more formal derivation of the instability?

Thank you in advance

Edit: I found this course about the finite volume method. But if there is something that explains it in the context of the finite element method or finite difference method it would be better for me.


问:
Finite difference/element method : time step and spatial resolution close to a finite singularity
I'm using the finite element method (FEM), but my question is quite a global question. It's related to this question but it is not the same.

Let's assume we have this equation : $$\partial_t c - u\frac{dc}{dx}=\frac{d}{dx}\Big(D(c)\frac{dc}{dx}\Big) $$

with $u$ a positive constant.

and with $D(c)>0$ and $$\lim_{c\rightarrow 1} D(c) \rightarrow{+\infty}$$

You can assume that there are Dirichlet conditions at the boundaries, for example : $c_{|\partial \Omega}=0$. And a Neumann boudary conditions : $c'(x)=0$ at the boundaries $\partial \Omega$.

What I'm observing is that for some reasons sometimes $c$ goes beyond the singularity $c=1$. And I don't want him to do that...

Here is a picture of the situation :

So I discretized the problem near the node $n$ where $c$ is close to 1 :

Let's first call $c_n^t$ the value of $c$ at node $n$ and time $t$.
The condition for being below $c=1$ is therefore : $\Delta t (u\frac{dc}{dx}+\frac{d}{dx}\Big(D(c)\frac{dc}{dx}\Big))<1-c_n^t$. If we discretize it, we get :

$$\Delta t(u\frac{c_{n+1}^t-c_{n-1}^t}{2\Delta x} + D(c_n^t)\frac{c_{n+1}^t+c_{n-1}^t-2c_{n}^t}{\Delta x^2 }+D'(c_n^t)(\frac{c_{n+1}^t-c_{n-1}^t}{2\Delta x})^2)<1-c_n^t$$

What I'm looking for is a general result like the one for the Péclet stability, where we need $Pe<1$. I'm talkinf about Péclet since here too we need $diffusion>convection$.

Thanks in advance


问:
Understanding MP-PIC implementation in OpenFOAM
The multiphase particle-in-cell (MP-PIC) method is characterized by mapping particle properties from the Lagrangian coordinates to the Eulerian grid. However, the implementation of this method in OpenFOAM solves an interpolated equation on the Eulerian grid only in one "submodel", the implicit packing model, for which there are alternatives (the none or explicit packing models). How is this possible and why is it still considered MP-PIC?

Furthermore, OpenFOAM seems to solve this equation in order to generate the particle velocity correction:
$$
\frac{\partial\alpha}{\partial t} - \Delta\left(\tau^\prime(\alpha) \alpha\right) = 0,
$$
where $\alpha$ is the particle phase volume fraction (averaged on the Eulerian grid) and $\tau$ is the interparticle stress. However, I can't find this equation in neither of the papers listed below. How does one derive the equation, and thus the velocity correction used in OpenFOAM?

References:

Andrews, M.J. and O'Rourke, P.J. (1996). The Multiphase Particle-in-Cell (MP-PIC) Method for Dense Particle Flows. International Journal of Multiphase Flow, 22(2):379–402.

Snider, D.M. (2001). An Incompressible Three-Dimensional Multiphase Particle-in-Cell Model for Dense Particle Flows. Journal of Computational Physics, 170:523–549.


问:
Interesting maxmin mixed integer/real quadratic optimization problem
I have the following problem:

$
\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}
\argmax_{\underset{\lambda_k\in \mathbb{R}}{\sigma_q^2(k)\in \mathbb{R}}}
\left[\min_{\underset{q\in \left[q_A,q_B\right]}{n_k\in \mathbb{Z}}}
\sum_{k=1}^N \frac{1}{\sigma_q^2(k)}\left(q-n_k\lambda_k\right)^2
\right] \qquad \mbox{s.t.}\qquad \sum_{k=1}^N  \frac{1}{\sigma_q(k)}=S
$

where the function to be minimized is quadratic both in the continuous variable $q$ and in the integer variables $n_k$ however I am having some problems translating it into a standard form. I tried to use the properties of the Hadamard product but that did not help much. Any ideas?


问:
Small residual but wrong results
When I use BiCGStab to solve a linear matrix system, I use the relative residual to exit the iteration and output the results. For calculating the relative residual I divide the norm of vector $r$ through the norm of vector $b$. Where $r$ is according to the definition of BiCGStab algorithm and $b$ is the right hand side of the matrix equation.

But even if I set the relative residual for convergency to 1e-8 the difference between calculated and exact results are quite big. When setting relative residual for convergency to 1e-14, I get good results.

The matrix systems I am solving have about 100k unknowns. A lot of entries in the result vector match the exact result. Only a few are off. But these few are enough to propagate the error.

Question: Is the relative residual really a good measure for convergency? Is this a normal behavior of the algorithm?

Title: Small residual but wrong results


问:
Monte Carlo domain not-so-dense
I already posted it on Physics SE, but maybe this is a better place:

I have a 5D integral being calculated with a Monte Carlo uniform random sampling. The issue is that the region of integration is very small and for 100000 points I get only around 20-30 points every time. The function has to be calculated several times for a certain parameter, which doesn't change the integration region but changes only the magnitude of the result (which is good, I suppose).

After several runs, I gathered some hundreds of points which are a valid match for the integration region. I was thinking about using those points to improve the integration result, but how can I do that, in a correct way?

EDIT: on the plane, the boundary is defined by an expression defining some kind of "little line" going from x = y = 150 up diagonally to infinity, so it's a thin domain. it is defined by all the five parameters in the function.

Title: Monte Carlo domain not-so-dense


问:
Implementing adaptive timestepping in CUDA
I want to implement a CUDA solver for the 2D shallow water equations using adaptive timestepping with a Courant number fixed by the user.  The algorithm pseudocode looks something like this:

while (t < end_time)
{
  evolve_flow_variables();
  dt = global_minimum_dt_over_all_elements_in_mesh();
  t += dt;
}
```

<p>I've seen example CUDA code with fixed timestepping (e.g. lecture 2 slides and code from Boston University) that makes multiple CUDA kernel calls from the time loops on the CPU, and synchronises only at the end:</p>

<pre class="lang-cpp prettyprint-override"><code>const dt;
while (t < end_time)
{
  evolve_flow_variables<<<grid_size, block_size>>>();
  t += dt;
}
cudaDeviceSynchronize();
</code>```

But I don't think this is feasible with adaptive timestepping.  Calling cudaDeviceSynchronize() after every timestep seems like a bad idea and, from my own measurements with Boston University's example code, the code could be up to twice as slow.

I'm aware of in-kernel grid synchronisation introduced in CUDA 9 that could allow me to put the entire time loop inside the kernel, but I'd rather be able to support slightly older devices.


问:
Is it possible to retain sets (physical tags) after remeshing the model in Gmsh?
I am currently  working on a remeshing algorithm for Abaqus using Gmsh as part of my Bachelors thesis and I seem to be stuck at some point.
I am able to remesh .geo and .msh files using the script given in Gmsh tutorial "t13".

t13.geo:

```<code>Merge "myJob5.msh";
angle = DefineNumber[40, Min 20, Max 120, Step 1,
  Name "Parameters/Angle for surface detection" ];
forceParametrizablePatches = 0;
forceIncludeBoundary = 0;
ClassifySurfaces{angle * Pi/180, forceIncludeBoundary, forceParametrizablePatches};
CreateGeometry;
Background Field = 1;
</code>```

The only problem I face is that my sets (physical tags) are still existent, but Gmsh didn't (or couldn't) match them with the appropriate nodes / lines / surfaces.
I started with a basic model.

myJob5.geo:

```<code>SetFactory("OpenCASCADE");
Point(1) = {0, 0, 0, 1.0};
Point(2) = {15, 0, 0, 1.0};
Point(3) = {15, 15, 0, 1.0};
Point(4) = {0, 15, 0, 1.0};
Point(5) = {7.5, 7.5, 0, 1.0};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line(5) = {1, 5};
Line(6) = {5, 3};
Line(7) = {2, 5};
Line(8) = {5, 4};
Physical Point("setNode1") = {3};
Physical Point("setNode2") = {2};
Physical Point("setNode2") += {1};
Physical Curve("setLine1") = {3};
Physical Curve("setLine2") = {4};
Physical Curve("setLine2") += {1};
Curve Loop(1) = {3, -8, 6};
Plane Surface(1) = {1};
Curve Loop(2) = {7, 6, -2};
Plane Surface(2) = {2};
Curve Loop(3) = {8, 4, 5};
Plane Surface(3) = {3};
Curve Loop(4) = {5, -7, -1};
Plane Surface(4) = {4};
Physical Surface("setSurf1") = {1};
Physical Surface("setSurf1") += {2};
Physical Surface("setSurf2") = {3};
Physical Surface("setSurf2") += {4};
</code>```

which translates into this picture:

After meshing, I ended up with this .msh file.

myJob5.msh:

```<code><span class="math-container">$MeshFormat
4.1 0 8
$</span>EndMeshFormat
<span class="math-container">$PhysicalNames
6
0 1 "setNode1"
0 2 "setNode2"
1 4 "setLine1"
1 5 "setLine2"
2 7 "setSurf1"
2 9 "setSurf2"
$</span>EndPhysicalNames
<span class="math-container">$Entities
5 8 4 0
1 0 0 0 1 2
2 15 0 0 1 2
3 15 15 0 1 1
4 0 15 0 0
5 7.5 7.5 0 0
1 -1.000000002804313e-07 -1e-07 -1e-07 15.0000001 1e-07 1e-07 1 5 2 1 -2
2 14.9999999 -1.000000002804313e-07 -1e-07 15.0000001 15.0000001 1e-07 0 2 2 -3
3 -1.000000002804313e-07 14.9999999 -1e-07 15.0000001 15.0000001 1e-07 1 4 2 3 -4
4 -1e-07 -1.000000002804313e-07 -1e-07 1e-07 15.0000001 1e-07 1 5 2 4 -1
5 -9.999999983634211e-08 -9.999999983634211e-08 -1e-07 7.5000001 7.5000001 1e-07 0 2 1 -5
6 7.4999999 7.4999999 -1e-07 15.0000001 15.0000001 1e-07 0 2 5 -3
7 7.4999999 -9.999999983634211e-08 -1e-07 15.0000001 7.5000001 1e-07 0 2 2 -5
8 -9.999999983634211e-08 7.4999999 -1e-07 7.5000001 15.0000001 1e-07 0 2 5 -4
1 -1.000000002804313e-07 7.4999999 -1e-07 15.0000001 15.0000001 1e-07 1 7 3 3 -8 6
2 7.4999999 -1.000000002804313e-07 -1e-07 15.0000001 15.0000001 1e-07 1 7 3 7 6 -2
3 -9.999999983634211e-08 -1.000000002804313e-07 -1e-07 7.5000001 15.0000001 1e-07 1 9 3 8 4 5
4 -1.000000002804313e-07 -9.999999983634211e-08 -1e-07 15.0000001 7.5000001 1e-07 1 9 3 5 -7 -1
$</span>EndEntities
<span class="math-container">$Nodes
12 5 1 5
0 1 0 1
1
0 0 0
0 2 0 1
2
15 0 0
0 3 0 1
3
15 15 0
0 4 0 1
4
0 15 0
0 5 0 1
5
7.5 7.5 0
1 1 0 0
1 3 0 0
1 4 0 0
2 1 0 0
2 2 0 0
2 3 0 0
2 4 0 0
$</span>EndNodes
<span class="math-container">$Elements
10 10 1 10
0 1 15 1
1 1
0 2 15 1
2 2
0 3 15 1
3 3
1 1 1 1
4 1 2
1 3 1 1
5 3 4
1 4 1 1
6 4 1
2 1 2 1
7 3 4 5
2 2 2 1
8 2 5 3
2 3 2 1
9 1 5 4
2 4 2 1
10 1 5 2
$</span>EndElements
</code>```

and this image:

Only the numbers are shown in this image to highlight the set that they are part of, e.g. the purple "0" means that that entity (surface, line or node) is not a part of any set, the yellow 4 means that the entity "line" is part of the set "setLine1".

Now, with the script of "t13" I tried to remesh the .msh file, giving me this image.

You can get this by simply 2d meshing the above given t13.geo file in gmsh.

Now my question is how I can retain the sets in my model, if possible at all with this method. What I am looking for is a way to match the newly generated elements with the previously assigned sets.

I know that this somehow must be possible since Gmsh does keep the sets accurate when refining the mesh, essentially multiplying the number of elements and assigning all of them to their respective set.

I really hope that my example is as clear and concise as possible. But still, if there is anything unclear, I'd love to clear things up and I am looking forward to anyone who can help me with this.
Thanks.


问:
Numerical Method for Equation System of two depending Equation Systems
I am searching a solution method for the following equation system of equation systems:

Let $A \in \mathbb{R}^{n \times n}$ be an invertible Matrix, $f, b_1, b_2 \in\mathbb{R}^n$ given vectors and $ C_i \in \mathbb{R}$ constants. I am searching for a numerical method to find the solution vectors $x_1, x_2  \in\mathbb{R}^n$ such that

\begin{align}
C_1A x_1 + C_2A x_2 &= C_1b_1 \\
C_3Ax_1 + C_4 Ax_2 &= f + C_3b_1  - C_5 b_2
\end{align}

I know that I also can write this as a equation system i.e.
\begin{equation}
\begin{pmatrix}
C_1 A & C_2 A \\
C_3 A & C_4 A
\end{pmatrix}
\begin{pmatrix}
x_1  \\
x_2
\end{pmatrix} =
\begin{pmatrix}
0 \\ f
\end{pmatrix} + \begin{pmatrix}
C_1 & 0\\
C_3 & -C_5
\end{pmatrix}\begin{pmatrix}
b_1 \\ b_2
\end{pmatrix}
\end{equation}

I am happy and looking forward getting some suggestions.


问:
Riemann solvers for metastable phases
Most Riemann solvers I've come across can solve the Riemann problem only under certain conditions such as convexity of the equation of state. But what happens if the fluid enters a metastable state or the spinodal region or the fluid state follows a nonconvex isentrope in the p-v diagram? All the approximate solvers described in the well known text books such as Toro break down here and I guess the same is true for exact solvers which only assume the three possible wave configurations of rarefaction,shock and contact.

But how to simulate for example a v.d.W - liquid including phase transitions?
What to do if the pressure of the unstable/metastable state becomes negative?
I've seen people doing the Maxwell construction in these cases. However, this results in in constant pressure-derivatives $\left(\frac{\partial p}{\partial \rho}\right)_S$ within these regions and as a result the adiabatic speed of sound becomes 0 which in turn makes the system of PDEs non-hyperbolic.

For the last couple of hours I was searching  for a more general Riemann solver. However, I wasn't able to find one.

Did I miss something or is this an unsolved problem?

Title: Riemann solvers for metastable phases


问:
BicgStab is not able to solve while Jacobi or GaussSeidel Methods can
I am trying to solve the 2D laplace equation,

$\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} = 0; \qquad   0 \lt x \lt 1, \quad0 \lt y \lt 1$

Subjected to the boundary conditions,

${\frac{d T}{d x}}|_{x=0,y}=0;$

${\frac{d T}{d y}}|_{x,y=0}=0;$

${\frac{d T}{d y}}|_{x,y=1}=0;$

$T|_{x=1,y} = 1; $

I am using Finite Difference method to discretise the differential equation and foward difference in the first two B.Cs and backward difference for the 3rd B.C. The solution which is $T=1$ at all grid points is perfectly coming out in Gauss Seidal as well as Jacobi iterative methods, but using BicgStab algorithm it is not converging even though the same matrices as Jacobi or Gauss Seidal. Any Ideas may help. Thank you.


问:
Book Recommendation: Analysis and design of mechanistic models - such as pharmacokinetics or hydrology models
I have been looking at an interesting book "Pharmacokinetic-Pharmacodynamic Modeling and Simulation" by Peter Bonate on pharmacokinetic models: the models of how medical drugs work their way through the tissues and organs of a human body, etc.

The book discusses "mechanistic" models where systems of differential equations model individual process in the body such as cellular metabolism or the drug, the operation of enzymes, etc. But these types of mechanistic models are used in many fields, such as hydrology to model flooding processes, weather forecasting, earthquake models, models for structural engineering to see how buildings will respond to different physical forces, etc.

While I have found some subject specific books on this topic, I was wondering if anyone knew of a good book or two on designing these types of mechanistic models? Given that these models have a relatively general form--such as large systems of equations--these models must face similar issues when it comes to calibrating against data, numerical methods, goodness of fit analysis, etc. dealing with multiple time-scales, etc.

Can anyone suggest a book that covers these topics across mechanistic models. Or even if someone know of good books within a particular subject area, that would be helpful. I was just hoping to standardize my own approaches to a problem I am working on.


问:
Finite element lemma proof
I was curious if anyone could help or provide a reference for the proof to the following lemma

Lemma:

Let $P_{1}$ be the set of polynomials of the first degree and let

$W = w(x) :  w \in C([0,1]), w(0) = w(1) = 0\;\text{and}\;w\mid_{(x_{i-1},x_{i})}$

then $W = V_{0}$ and

$w(x) = \sum_{i=1}^{n-1} w(x_{i}) \phi_{i}(x) \forall w(x) \in W = V_{0}$
To note, $V_{0} = \text{span}(\phi_{1},\dots,\phi_{n})$

Title: Finite element lemma proof


问:
Ramp least squares estimation
With some given $s$ value, let
\begin{equation}
\begin{aligned}
h(\beta)&=\min(\sum_{i=1}^n(Y_i - X_i\beta)^2, s)\\
        &=\sum_{i=1}^n(Y_i - X_i\beta)^2-\max(0, \sum_{i=1}^n(Y_i - X_i\beta)^2-s)\\
        &=f(\beta)-g(\beta, s)
\end{aligned}
\end{equation}
The problem is:
\begin{equation}
\min h(\beta)
\end{equation}
I used Disciplined Convex-Concave Programming and DCCP package based on CVXPY. The code follows:

import cvxpy as cp
import numpy as np
import dccp
from dccp.problem import is_dccp

# Generate data.
m = 20
n = 15
np.random.seed(1)
X = np.random.randn(m, n)
Y = np.random.randn(m)

# Define and solve the DCCP problem.
def loss_fn(X, Y, beta):
    return cp.norm2(cp.matmul(X, beta) - Y)**2

def obj_g(X, Y, beta, sval):
    return cp.pos(loss_fn(X, Y, beta) - sval)

beta = cp.Variable(n)
s = 10000000000000
constr = obj_g(X, Y, beta, s)
t = cp.Variable(1)
t.value = [1]
cost = loss_fn(X, Y, beta) - t
problem = cp.Problem(cp.Minimize(cost), [constr >= t])
print("problem is DCP:", problem.is_dcp())   # false
print("problem is DCCP:", is_dccp(problem))  # true

problem.solve(verbose=True, solver=cp.ECOS, method='dccp')
# Print result.
print("\nThe optimal value is", problem.value)
print("The optimal beta is")
print(beta.value)
print("The norm of the residual is ", cp.norm(X*beta - Y, p=2).value)
```

<p>Because of the large value $s$, I would hope to get a solution similar to the least squares estimation. But there is no solution as the output shows:</p>

<pre class="lang-py prettyprint-override"><code>problem is DCP: False
problem is DCCP: True

ECOS 2.0.7 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS

It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
 0  +0.000e+00  -0.000e+00  +2e+01  9e-02  1e-04  1e+00  9e+00    ---    ---    1  1  - |  -  -
 1  -7.422e-04  +2.695e-09  +2e-01  1e-03  1e-06  1e-02  9e-02  0.9890  1e-04   2  1  1 |  0  0
 2  -1.638e-05  +5.963e-11  +2e-03  1e-05  2e-08  1e-04  1e-03  0.9890  1e-04   2  1  1 |  0  0
 3  -2.711e-07  +9.888e-13  +2e-05  1e-07  2e-10  2e-06  1e-05  0.9890  1e-04   4  1  1 |  0  0
 4  -3.991e-09  +1.379e-14  +2e-07  1e-09  2e-12  2e-08  1e-07  0.9890  1e-04   1  0  0 |  0  0
 5  -5.507e-11  +1.872e-16  +3e-09  2e-11  2e-14  2e-10  1e-09  0.9890  1e-04   1  0  0 |  0  0

OPTIMAL (within feastol=1.6e-11, reltol=4.8e+01, abstol=2.6e-09).
Runtime: 0.001112 seconds.

ECOS 2.0.7 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS

It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
 0  +0.000e+00  -5.811e-01  +1e+01  6e-01  6e-01  1e+00  2e+00    ---    ---    1  1  - |  -  -
 1  -7.758e+00  -2.575e+00  +1e+00  2e-01  7e-01  6e+00  3e-01  0.9890  1e-01   1  1  1 |  0  0
 2  -3.104e+02  -9.419e+01  +4e-02  2e-01  8e-01  2e+02  8e-03  0.9725  8e-04   2  1  1 |  0  0
 3  -2.409e+03  -9.556e+02  +5e-03  2e-01  8e-01  1e+03  1e-03  0.8968  5e-02   3  2  2 |  0  0
 4  -1.103e+04  -5.209e+03  +2e-03  2e-01  7e-01  6e+03  4e-04  0.9347  3e-01   2  2  2 |  0  0
 5  -1.268e+04  -1.592e+03  +8e-04  1e-01  1e+00  1e+04  2e-04  0.7916  4e-01   3  2  2 |  0  0
 6  -1.236e+05  -2.099e+04  +9e-05  1e-01  1e+00  1e+05  2e-05  0.8979  9e-03   1  1  1 |  0  0
 7  -4.261e+05  -1.850e+05  +4e-05  2e-01  7e-01  2e+05  1e-05  0.7182  3e-01   2  1  1 |  0  0
 8  -2.492e+07  -1.078e+07  +7e-07  1e-01  7e-01  1e+07  2e-07  0.9838  1e-04   3  2  2 |  0  0
 9  -2.226e+08  -9.836e+07  +5e-08  9e-02  5e-01  1e+08  1e-08  0.9339  2e-03   2  3  2 |  0  0

UNBOUNDED (within feastol=1.0e-09).
Runtime: 0.001949 seconds.

The optimal value is None
The optimal beta is
None
The norm of the residual is  None
</code>```

Title: Ramp least squares estimation


问:
Eigenfaces Algorithm
This might be a silly quesntion but recently I've been trying to program the eigenface algorithm using PCA, so I arranged the face vectors vertically in a matrix X such as: X = [x1,x2,x3,...,xn];
In this case, what would be the wright way of computing the covariance matrix? cov(X) or cov(X')?</p>

Title: Eigenfaces Algorithm


问:
Palmprint Identification - Why do we align the images before we use the Fourier Transform?
I am currently reading the paper PALMPRINT IDENTIFICATION BY FOURIER TRANSFORM by WENXIN LI, DAVID ZHANG and ZHUOQUN XU about identifying persons based on an image of their palm, one version of the paper can be found here.

The first step in their identification algorithm consists of aligning and preprocessing the image - since they spent a few pages on that front, this gives me the impression, that this step is kinda important.

What I don't understand however: Why do we do this? The essential part in my opinion is the Fourier Transform afterwards, that is used to extract/generate the feature. But from my understanding, this should even work without this aligning/preprocessing step.

I have found a property of the DFT that states that rotating the spatial domain contents rotates the frequency domain contents, can be read here.
In order to extract features, some rotation might be necessary, at some point. In my opinion however, this does not justify to think up such an elaborate method of aligning/preprocessing the image, there has to be more behind that.

Therefore, I am most likely overlooking/missing something. Thanks in advance.


问:
Partial/Extended/Truncated Template Matching
So template matching using correlation is available in a lot of computational packages; OpenCV matchTemplate(), scipy.signal.correlate2d(), IPP CrossCorrNorm, etc.

But they all either evaluate offsets for which the template is fully inside the image, or use a boundary condition like padding or wrapping to evaluate offsets where the template hangs outside the image. I don't like that, because in general there's no good reason to believe those various kinds of pads are going to correlate with the overhanging portion of the template.

Is there any implementation out there which includes a boundary condition like 'truncate', where normal 'full' correlation is computed for all template offsets within the image, but also for all offsets where the template sticks out, the correlation is computed only for the subset of template pixels within the image?

At the extreme offsets, the template and image would share only one pixel, at which point normalized correlation is undefined (because of divide by standard deviation of 1 sample = 0), but those corners of the correlation surface can just be 0.

The benefit would be to be able to match a template even if it appears only partially in the image (see illustration below). It would be up to the caller to determine how close to the edge of the extended correlation surface (what fraction of the template) to accept as a match.

For an Image of size HxW and a template of size hXw, the resulting full correlation surface is (H-h+1)x(W-w+1), but for this extended/dynamic/partial template matching the correlation surface would be of size (H+h-1)x(W+w-1) (and again, the corners would be 0 because they are undefined)

I have an implementation of this, it's not that hard, but it's also not efficient. Well-optimized image correlation takes good advantage of repeated sums/products. Before I spend effort on optimizing my own implementation, it would be good to know if there is already a good implementation out there.


问:
What algorithm do BLAS and ATLAS use for matrix multiplication?
I have searched and what I understood was that they use the naive one with several memory and cache optimizations.

However, I wanted to know whether they are using the Strassen or the Coppersmith-Winograd algorithms.

If they don't use them, then why?

Additionally, does NumPy use BLAS or ATLAS? What I did is that I looked up similar questions but did not find specific answer to why/why not NumPy is faster than most implementations.


问:
implementation for coppersmith matrix multiplication
Is there any online implementation for the coppersmith matrix multiplication I have searched alot but can not find any?  and if there is not any why is that Isn't this algotithm much faster than strassen in very large square matrices in the range of 100000 to millions,  from what I understand the only downside of it is that  the constants of running time is high</p>


问:
Smoothing FFT result
I am trying to calculate the spectrum of Bremmstrahlung, which involves calculating the Fourier transformed acceleration. I am solving a non-linear ODE to numerically calculate the acceleration in the time domain. After taking the Fourier transformation using Numpy's fft, the resultant spectrum looks highly non-smooth and "non-physical" .
I cannot paste the entire code so I am posting what I think is relevant snippet. Can someone point out what I am doing wrong?

Note: My acceleration is a function of two variables (beta, and b, the impact parameter), and I want to plot it for the different b, and in order to factor out the beta term I am just summing over all the values of acceleration for the different beta, for a given impact parameter b.

Also my spectrum is Fourier transformed acceleration square times a constant factor (Larmor's formula)

```<code>'''
Fourier Transformation of the acceleration
'''
N = 2**8
sampling_frequency = 10
a_w_normalized = [[] for a in range(len(impact_parameter))]
a_w =[[] for a in range(len(impact_parameter))]
intensity_normalized = [[] for a in range(len(impact_parameter))]
acc_summed_over_velocity = []
intensity_summed_over_velocity =[]

for index,b in enumerate(impact_parameter):
    acc_sum =[0 for x in range(N)]
    intensity_sum = [0 for x in range(N)]
    for j in range(len(velocity_z_component)):
        #window_kaiser = signal.kaiser(N, 15)
        #window_hann = signal.hann(N,sym=True)
        window =1
        fft_input = acc_normalized[index][j]*window
        ft_acc_normalized = np.abs(np.fft.fft(fft_input,norm=None))

        acc_sum =np.add(ft_acc_normalized,acc_sum )
        intensity_list = [power_spectrum_factor * (a ** 2) for a in ft_acc_normalized]
        intensity_sum = np.add(intensity_list,intensity_sum)
        intensity_normalized[index].append(intensity_list)
        a_w_normalized[index].append(ft_acc_normalized)
    acc_summed_over_velocity.append(acc_sum)
    intensity_summed_over_velocity.append(intensity_sum *acceleration_factor**2)

#intensity_summed_over_velocity=  ma.masked_less_equal(intensity_summed_over_velocity,1e-5)

'''
Plotting acceleration for selected values of impact paramters in time and frequency domain
'''

if plot is True:
    plt.figure(figsize=(12, 8))
    for i in range(len(impact_parameter)):
        plt.plot(t,acc_timedomain_summed_over_velocity[i], label='b={:.3f}'.format(impact_parameter[i]), )
    plt.legend()
    plt.ylabel(r'<span class="math-container">$ a(\tilde t) $</span>', fontsize=14)
    plt.xlabel(r'<span class="math-container">$ \tilde t $</span>', fontsize=14)
    ax.spines['left'].set_position('center')
    ax.spines['bottom'].set_position('center')
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['left'].set_smart_bounds(True)
    ax.spines['bottom'].set_smart_bounds(True)
    ax.xaxis.set_ticks_position('bottom')
    ax.yaxis.set_ticks_position('left')
    plt.xticks(fontsize=15)
    plt.yticks(fontsize=15)
    if screening is False:
        plt.title('a(t) vs time without screening',fontsize=15)
        plt.savefig('./Plots/acc_without_screening.png')
    else:
        plt.title('a(t) vs time with screening',fontsize=15)
        plt.savefig('./Plots/acc_with_screening.png')

    '''
    Fourier transformed intensity for different impact paramater
    '''

    plt.figure(figsize=(12, 8))
    freq_normalized = np.fft.fftfreq(N)*(2*np.pi*sampling_frequency)
    for index,b in enumerate(plasma.impact_parameter):
        plt.plot(np.abs(freq_normalized), intensity_summed_over_velocity[index], label=r'<span class="math-container">$ \tilde b={:.2f}$</span>'.format(b), )
    plt.xlabel(r'<span class="math-container">$ \tilde \omega $</span>', fontsize=14)
    plt.ylabel(r'<span class="math-container">$ I_{\omega} $</span>', fontsize=16)
    plt.legend(loc="lower left")
    plt.xscale('log')
    plt.yscale('log')
    plt.xticks(fontsize=12)
    plt.yticks(fontsize=12)

    if screening is False:
        plt.title('Single particle spectrum without screening')
        plt.savefig('./Plots/spectrum_no_screening.png')
    else:
        plt.title('Single particle spectrum with screening')
        plt.savefig('./Plots/spectrum_debye_screening.png')

    plt.show()
</code>```

Title: Smoothing FFT result


问:
Forming Basis Functions from 6-31G Basis Set for Carbon Atom
I am a computer science grad and I am working to write an electronic structure calculation program and I am stuck at forming basis functions using 6-31G Basis set for atoms having higher atomic numbers (p,d,f shells).

For example, consider a gaussian g as

```<code>g(alpha) = (2*alpha/pi)^(3/4) * exp(-alpha*r^2).
</code>```

Looking at the Carbon (1s2 2s2 2p2) basis set data for 6-31G basis set

```<code>C     0
S   6   1.00
      0.3047524880D+04       0.1834737132D-02
      0.4573695180D+03       0.1403732281D-01
      0.1039486850D+03       0.6884262226D-01
      0.2921015530D+02       0.2321844432D+00
      0.9286662960D+01       0.4679413484D+00
      0.3163926960D+01       0.3623119853D+00
SP   3   1.00
      0.7868272350D+01      -0.1193324198D+00       0.6899906659D-01
      0.1881288540D+01      -0.1608541517D+00       0.3164239610D+00
      0.5442492580D+00       0.1143456438D+01       0.7443082909D+00
SP   1   1.00
      0.1687144782D+00       0.1000000000D+01       0.1000000000D+01
</code>```

How many basis functions will carbon have in this case? I have been following a lecture series which does HF SCF calculation for a simple H2 molecule where each Hydrogen atom has two basis functions as per 6-31G basis set.

```<code>H     0
S   3   1.00
      0.1873113696D+02       0.3349460434D-01
      0.2825394365D+01       0.2347269535D+00
      0.6401216923D+00       0.8137573261D+00
S   1   1.00
      0.1612777588D+00       1.0000000
</code>```

I can write the two basis functions in the form

```<code>b1= 0.0334*g(18.73) + 0.2347*g(2.825) + 0.8137*g(0.64);
b2= g(0.161);
</code>```

In a similar manner, how would I go about writing the basis functions for Carbon atom. Any advise would be very helpful as I am stuck in this for couple of days and haven't been able to find any substantial resource for this.

I did stumble upon this previous post but having a hard time figuring out how did 9 basis functions come for Carbon.


问:
Setting up diffusion with integral B.C. in Fenics
I'm trying to model diffusion through a cylindrical domain $D = \{ (x,y,z) : x^2 + y^2 \leq 1, \;\; 0 \leq z \leq 1\}$.

The is an initial concentration of the diffusant at the upper flat surface, say $C(x,y,z=1) = C_{init}$, and initially zero concentration within the domain.

The diffusant can move through both the $z=0$ and $z=1$ interfaces; there is no flux through the vertical walls.

The diffusion equation is
$$ \frac{\partial C}{\partial t} = \Delta C. $$

The boundary conditions are the problem:
$$ \frac{\partial C}{\partial t} (x, y, z=0, t) = K_0 \int_{z=0 \text{ cap}} \frac{\partial C}{\partial z} \;dS,
$$
$$ \frac{\partial C}{\partial t} (x, y, z=1, t) = K_1 \int_{z=1 \text{ cap}} \frac{\partial C}{\partial z} \;dS.
$$
These come from assuming that the diffusion through the end caps occurs into finite-volume chambers, and the accumulated diffusant affects the rate of diffusion across the boundary. The $K_i$ are constants.

Trying to set up the variational formulation for Fenics, with $C_N$ equal to the concentration at time step $N$, I get

$$ \int_D \frac{C_{N} - C_{N-1}}{\Delta t} v \; dV + \int_D \nabla C_N \cdot \nabla v \; d V - \int_{\partial D} v \nabla C_N\; d S = 0, $$

or equivalently
$$ \int_D \frac{C_{N} - C_{N-1}}{\Delta t} v \; dV + \int_D \nabla C_N \cdot \nabla v \; d V - \int_{z=1 \text{ cap}} v  \frac{\partial C_N}{\partial z}\; d S + \int_{z=0 \text{ cap}} v  \frac{\partial C_N}{\partial z}\; d S= 0. $$

I don't know where to go from here. The boundary conditions tell about the integral of $\partial C / \partial z$ over the caps, but I don't know how to translate that into information about the integral against $v$ over the caps. Any advice would be welcome.

N.B.: I know that the problem I've stated is essentially one dimensional, and can be solved that way. The problem I'm actually interested in has a more complex domain, so I gave a simpler version that should capture the B.C. issue.


问:
Get interpolated values in user defined 2D grid Paraview
I have a 2D flow and would like to obtain the value of certain scalar field in a set of points forming a regular mesh. These points should not coincide with the nodes of the actual mesh used in the CFD simulation.
May be you know how to automatize this process in ParaView.

Basicly, I need to first define the set of points, get the value of u at them, and save them in a matrix with x and y direction.

Cheers


问:
Is there an open-source material database management GUI?
Does somebody know an open-source GUI for the management of a small material database?

I have a spreadsheet with some materials in it. Each materials has some temperature-dependent properties like this:

.

The problem is that some parameters are a function of temperature, and entering equations in cells is a little error-prone. Also I would like to switch between units (imperial, SI) sometimes.

I thought about writing a small GUI for this (perhaps in Python or a Jupiter notebook, ideas welcome :) ), but before I reinvent the wheel, I wanted to ask you guys if there is already something like this.

Cheers,

Ken


问:
Determining the indices of a VTK mesh
I am trying to assign fixed constraint at specific indices of VTK mesh, however, I only can view STL file in a blender as follows:

Upon assigning the favored indices in my scene at the SOFA physics simulator, of course, I did not get favored results since I am using VTK mesh.

To view VTK mesh, I am using GMSH, the mesh sounds cloudy and I don't know what could be the easiest way to identify the vertices, I aiming to fix one end of the twisted beam simply

Thanks in advance

Title: Determining the indices of a VTK mesh


问:
How to ensure values stay within range?
e.g. water in a height map

Choosing a range with a margin of error for typical model behaviour seems practical.
Could we instead (1). predict maximum values;  or (2). have a natural maximum?

1. predict maximum The total energy in a system probably limits the maximum possible values (i.e. if it was all concentrated into one single value).
For example, a ring of waves moving towards a center (opposite of ripples on a puddle) would superpose at the center, summing their heights.

Or, for a lower limit, could we calculate the directions of waves and their energy, and find the maximum possible collisions (or are there no known tricks to shortcut its computational complexity)? Or even, predict the actual collisions?

2. natural maximum The speed of light is a natural limit; as you approach it, the same force produces less acceleration. Relativity is a consistent system for nature - so relativistic fluid dynamics should be too. We could choose a different maximum speed, perhaps walking pace. This doesn't seem useful for modelling real systems, but it would limit the range of values in a "natural" way.

(Perhaps $c$ is a way to avoid instability in the simulation we inhabit).

Title: How to ensure values stay within range?


问:
Solving Parabolic PDE using Matlab
I have the following pde (Burger's equation)

for $\epsilon>0: u_t+u.u_x=\epsilon.u_{xx}$ and $x\in \mathbb{R},t>0$

and the initial condition: $u(x,0)=\phi(x)=1_{(-\infty,0)}(x)$.

I want first to solve it (using matlab) for the case $\epsilon=0$, which is a first pde, and then for $\epsilon=0.1$ and $ \epsilon=0.001$.

Title: Solving Parabolic PDE using Matlab


问:
Extracting FEM matrices in matlab pde toolbox
I am trying to follow the dynamic linear elasticity in Matlab, link here. My goal is to extract the FE Matrices using the function assembleFEMatrices in matlab and solve the resulting system of second-order ODE's by Backward Euler, for example.

I am not having any success reconciling the result from solvepde with the result from my code.

Let ''model'' be the PDE model specified according to the Matlab example.

To extract the FE Matrices, I do:

```<code>FEM = assembleFEMatrices(model);
M = FEM.M;
K = FEM.K;
F = FEM.F;
G = FEM.G;
</code>```

This is following Matlab's documentation on FEM, on writing the PDE, and on extracting the FE Matrices.

Let $x,y$ be the displacement in the $x,y$ directions, respectively. The resulting system of ODE's should be

$M \begin{bmatrix} \ddot{x} \\ \ddot{y} \end{bmatrix} + K \begin{bmatrix} x \\ y \end{bmatrix} = G+F$.

Here, $M$ is the mass matrix, $K$ is the stiffness matrix, $G$ is from the Neumann BC while $F$ is for the forcing term.

I then write the above system as

$ \begin{bmatrix} \ddot{x} \\ \ddot{y} \end{bmatrix} + M^{-1}K \begin{bmatrix} x \\ y \end{bmatrix} = M^{-1}(G+F)$.

Then to reduce this to a 1st-order ODE, I introduce $z_1 = \begin{bmatrix} x \\ y \end{bmatrix}$ and $z_2 = \begin{bmatrix} \dot{x} \\ \dot{y} \end{bmatrix}$

to obtain

$\frac{d}{dt} \begin{bmatrix} z_1 \\ z_2 \end{bmatrix} = \begin{bmatrix} 0 & I \\ -M^{-1}K & 0 \end{bmatrix} \begin{bmatrix} z_1 \\ z_2 \end{bmatrix} + \begin{bmatrix} 0 \\ M^{-1}(G+F) \end{bmatrix}$.

I then solve the above system using backward euler and compared $z_1$ with the result I get from Matlab's solvepde but they don't match. I don't think the problem is with my solution approach because I tried Matlab's lsim and I get almost the same results as backward euler. Did I incorrectly extract and misinterpret the FE Matrices?

Suggestions appreciated.

Edit: Code below. I slightly changed some constants in the Matlab example for the code to run faster. I also considered the full beam instead of just half the beam. In my code below, mySoln is the solution I get from assembleFEMatrices and transforming the 2nd order ODE to a first order ODE as I described above. results.NodalSolution is the solution obtained from Matlab's solvepde. I can't get this to reconcile with mySoln.

```<code>% Create PDE model

N = 2; % Two PDE in plane stress elasticity
model = createpde(N);
blength = 1.5; % Beam length, in.
height = .01; % Thickness of the beam, in.

l2 = blength/2;
h2 = height/2;
rect = [3 4 0 l2 l2 0 -h2 -h2  h2 h2]';
g = decsg(rect,'R1',('R1')');
pg = geometryFromEdges(model,g);

figure
pdegplot(g,'EdgeLabels','on');
title('Geometry With Edge Labels Displayed');
axis([-.1 1.1*l2 -5*h2 5*h2]);

E = 3.0e7; % Young's modulus of the material, lbs/in^2
gnu = .3; % Poisson's ratio of the material
rho = .3/386; % Density of the material
G = E/(2.*(1+gnu));
mu = 2*G*gnu/(1-gnu);
c = [2*G+mu; 0; G;   0; G; mu; 0;  G; 0; 2*G+mu];
f = [0 0]'; % No body forces
specifyCoefficients(model, 'm', rho, 'd', 0, 'c', c, 'a', 0, 'f', f);

% BCs
symBC = applyBoundaryCondition(model,'dirichlet','Edge',4,'u',[0 0]);
clampedBC = applyBoundaryCondition(model,'dirichlet','Edge',2,'u',[0 0]);
sigma = 5e1;
presBC = applyBoundaryCondition(model,'neumann','Edge',3,'g',[0 sigma]);
setInitialConditions(model,0,0);
generateMesh(model,'Hmax',height/2,'MesherVersion','R2013a');

% Solve pde using matlab pde toolbox
tfinal = 3e-3;
tlist = linspace(0,tfinal,100);
result = solvepde(model,tlist);

%%
% plot appearance of beam

figure;

X_nodes = result.Mesh.Nodes(1,:);
X_nodes = X_nodes(:);
Y_nodes = result.Mesh.Nodes(2,:);
Y_nodes = Y_nodes(:);

for k = 1:size(result.SolutionTimes,2)

    % Extract displacement at time k
    disp_X = result.NodalSolution(:,1,k);
    disp_X = disp_X(:);

    disp_Y = result.NodalSolution(:,2,k);
    disp_Y = disp_Y(:);

    plot(X_nodes + disp_X, Y_nodes + disp_Y,'*')
    pause(0.1)

end

%%

% try to recover solution by assembling the FEM matrix and time-stepping it
FEM = assembleFEMatrices(model);
Mmat = FEM.M;
Kmat = FEM.K;
Fmat = FEM.F;
Gmat = FEM.G;

nNodes = size(model.Mesh.Nodes,2);

Amat = [zeros(2*nNodes,2*nNodes) eye(2*nNodes); -Mmat\Kmat zeros(2*nNodes,2*nNodes)];
Bmat = [zeros(2*nNodes,1); Mmat\(Fmat + Gmat)];

% Do backward euler here

deltaT = tlist(2) - tlist(1);
IminushA = sparse(eye(4*nNodes)) - deltaT*Amat;

mySoln = zeros(4*nNodes,length(tlist));

for k = 2:length(tlist)
   fprintf('%d...\n',k)
   mySoln(:,k) = IminushA\(mySoln(:,k-1) + deltaT*Bmat);
end
</code>```


问:
Solving diffusion equation using finite difference method
I am solving an 1-dimensional diffusion equation with Neumann boundary condition at outlet and constant concentration, C, at the inlet. In the end, I want to observe how the concentration diffuses over time along the x direction.

$$\frac{\partial C}{\partial t} =  D\frac{\partial ^2 C}{\partial x^2}$$

The following is the code written in MATLAB.

```<code>function main()
C0 = [500 480 460 440 420 400 380 360 340];
tspan = [0 10];
[t C]  = ode23s(@(t,s) fun(t,s), tspan , C0);
plot(t,C)

function dC = fun(t,C)
D  = 3000;
retain_rows = 2:8;
tri = -(D/(1/5^2))*full(gallery('tridiag',9,-1,2,-1)); % 5 is del x
tri = tri(retain_rows,:);
dC(1,1) = 0; % constant concentration at inlet
dC(retain_rows,1) = tri*C;
dC(9,1) = -D*((C(9)-C(8))/2); Neumann boundary condition at outlet
end
end
</code>```

I have the solution obtained from ode solver. I'd like to know how to check for mole balance from the solution that is obtained.
Also, could someone explain how the above code can be changed to apply constant flux at the inlet?

Any suggestions?


问:
How to use matlab **fft** function or other fast methods to solve a convection-diffusion system quickly?
for a convection-diffusion equation with Dirichlet boundary conditions as follows:
$$-u''+qu'=f$$
Using centered difference for $u''$ and $u'$, we get a linear system
$$Ax=b$$where matrix $A$ is nonsymetric.
Let $H = (A+A')/2,S = (A-A')/2$, we know that matrix $H$ is symmetric positive definite. Given a constant $\alpha>0$, we can get two linear large sparse systems, and $I$ is a identity matrix.
$$(\alpha I+H)x_{k+1/2} = b_1$$and
$$(\alpha I+S)x_{k+1} = b_2$$

My question is in 3D case, the above two systems are so large and sparse, how to use sin and modified sin transfomation to solve the 2 systems quickly, or some other fast methods, because in matlab, the backslash command '\' tend to be slow for this very large sparse in 3D. Can someone give the matlab codes, thanks very much.

Below is my matlab codes: for n=64 equidistance so slow.

clc;clear;close all
n=64;q=1;
h=1/(n+1);
I = speye(n);
r = q*h/2;
t1 = 6; t2 = -1-r; t3 = -1+r;
Tx = spdiags(ones(n,1).*[t2 t1 t3],[-1 0 1],n,n);
Ty = spdiags(ones(n,1).*[t2 0 t3],[-1 0 1],n,n);
Tz = spdiags(ones(n,1).*[t2 0 t3],[-1 0 1],n,n);
A = kron(Tx,kron(I,I))+kron(I,kron(Ty,I))+kron(I,kron(I,Tz));
H = (A+A')/2;   S = (A-A')/2;
lambda_maxH = eigs(H,1,'largestabs');
lambda_minH = eigs(H,1,'smallestabs');
alpha = sqrt(lambda_maxH*lambda_minH);
N = length(A);
I = speye(N);
b = A*ones(N,1);

%   this 2 systems are solved so slow
x1 = (alpha*I+H)\b;
x2 = (alpha*I+S)\b;

```


问:
Efficient evaluation of spherical harmonic expansions
Assume I know that I can express an approximation of a function by
$$
\sum_{l=0}^{k}\left( \sqrt{A_l} z_{l,0}^1 L_{l,0}(\theta)+\sqrt{2A_l}\sum_{m=1}^l L_{l,m}(\theta)(z_{l,m}^1 \cos(m\phi)+z_{l,m}^2\sin(m\phi))\right)
$$
where $\theta$ is the inclination and $\phi$ azimuth, and $z_{l,m}^{1,2}$ are a sequence of indepedent standard normal random numbers, and $A_l$ is some pre-decided number so that it decays fast enough. The functions $L_{l,m}$ are given by
$$
L_{l,m}=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}P_{l,m}(\cos(\theta))
$$
where $P_{l,m}$ are the associated Lagrange polynomial

of the first type.  This I know. I would like to evaluate this sum for

1) Large values of the truncation parameter $k$ (like 10^6)

2) many many values of $\theta$ and $\phi$.

I would like to do this in python, but i suppose I would be ok with writing a compiled program in Fortran for instance, as long as I can acess from python.

My attempt so far is to use scipy.special.lpmv and just blindly evaluate the sum, but that breaks down very fast for not very large values of l,m.
There seems to be some sort of trouble with the scipy implementation, see.

https://stackoverflow.com/questions/39249594/numerical-computation-of-the-spherical-harmonics-expansion

```<code>import numpy as np
from scipy.special import lpmv sph_harm
from math import factorial
from dolfin import *
from mshr import *
from numpy import sqrt,cos,sin,pi
import numpy as np
no_terms=50
def rand_sequence(noterms,seed):
    np.random.seed(seed)
    return(np.random.normal(0,1,[noterms,noterms,2]))
def angular_spec(l,alpha):
    return((1+l)**(-alpha))
def plm(l,m,theta):
    return(sqrt((2*l+1)*factorial(l-m)/4/pi/factorial(l+m))*lpmv(m,l,cos(theta)))    #it seems to break here!!
rands=rand_sequence(no_terms,123)

def GRF(noterms,alpha,rands,x,y,z):
    phi=np.arctan2(y,x)
    theta=np.arccos(z)
    out=0.0
    for l in range(noterms):
        a=sqrt(angular_spec(l,alpha))*rands[l,0,0]*plm(l,0,theta)
        b=sqrt(2*angular_spec(l,alpha))*sum([plm(l,m,theta)(rands[l,m,0]*cos(m*phi)+rands[l,m,1]*sin(m*phi)) for m in range(1,l+1)])
        out=np.add(out,np.add(a,b))
    return(out)
</code>```

So, any suggestions? Will I have to give up my hope of being able to evaluate for large $k$? How can I speed this up. I will have a very large amount of $\theta$s and $\phi$s to evaluate in...


问:
Live audio processing using c++ (standard or external libs 'ue4')
My goal is to get audio from the input of another device (either connected to computer throught aux or audio interface threw rca or coaxil(from an external device playing audio live) adapters to line or something) and let a certain property (gain or volume, probably just volume) drive the animation of a 3d model (The more the volume increases the model moves along a directional axis 'either procedraly animated or let it drive a pre made animation) preferably all in unreal engine because it already has a realtime graphics engine so no need to worry about the graphics just how would I go about getting the signal (which libraries should I look into) and how would I let the data drive the animation (not specific code that would work specifically in unreal just broad methodology)

Thanks a bunch


问:
Hybrid Ellpack-Itpack (ELL) + COO Sparse Matrix Representation decomposition threshold
Hybrid ELL-COO sparse matrix representation can be done as in the picture, I was looking intensively, however I couldn't find out what is the threshold of decomposing the original matrix into ELL part and COO part, how to divide the original matrix A ?


问:
Probability approximation: monte carlo VS sde
I have a probability measure $\mu$ (say, in $\mathbb{R}^{d}$, with density) and I want to approximate it numerically. Today I noticed that my measure is ergotic for a certain Stochastic Differential Equation. If we call $x$ our starting point (asymptotically not relevant), to sample from $\mu$ I see two alternatives:

run the usual Metropolis-Hasting chain (MCMC);

simulate a SDE path starting from $x$ (e.g. via Euler's discretization), then stop at large times;

I am  right, or did I misunderstand the strategy number 2?
How is the efficiency of these two algorithms compared, in general?

My book on Monte Carlo methods does not consider SDEs, while mine on SDEs does ignore the Metropolis algorithm. Good rigorous references (as well as papers) are highly appreciated. Thanks in advance.


问:
How to handle system of chemical reactions for a batch reactor SciPy solver
I have a system of chemical reactions where the rate equations represent a batch reactor model. The model is a system of ODEs which is solved with the SciPy solve_ivp function. My example (see below) currently contains 5 reactions but I need to add many more reactions from the overall scheme. With my current approach it has been tedious to keep track of the chemical species, mass fractions, rate parameters, etc.

```<code>import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Parameters
# -----------------------------------------------------------------------------

species = (
    'CELL', 'CELLA', 'CH2OHCHO', 'CHOCHO', 'CH3CHO', 'C6H6O3', 'C2H5CHO',
    'CH3OH', 'CH2O', 'CO', 'GCO', 'CO2', 'H2', 'H2O', 'GCH2O', 'HCOOH',
    'CH2OHCH2CHO', 'CH4', 'GH2', 'CHAR', 'C6H10O5', 'GCOH2', 'GMSW', 'HCE1',
    'HCE2'
)

# time vector to evaluate reaction rates [s]
time = np.linspace(0, 2.0)

# Solve batch reactor system of equations
# -----------------------------------------------------------------------------

def dcdt_debiagi(t, y):
    """
    Rate equations for Debiagi 2018 biomass pyrolysis kinetics.
    """
    R = 1.9859          # universal gas constant [cal/(mol K)]
    T = 773.15          # temperature [K]

    cell = y[0]         # cellulose mass fraction [-]
    cella = y[1]        # active cellulose mass fraction [-]
    gmsw = y[22]        # softwood hemicellulose mass fraction [-]

    # Cellulose reactions and rate constants
    # 1) CELL -> CELLA
    # 2) CELLA -> 0.40 CH2OHCHO + 0.03 CHOCHO + 0.17 CH3CHO + 0.25 C6H6O3 + 0.35 C2H5CHO + 0.20 CH3OH + 0.15 CH2O + 0.49 CO + 0.05 GCO + 0.43 CO2 + 0.13 H2 + 0.93 H2O + 0.05 GCH2O + 0.02 HCOOH + 0.05 CH2OHCH2CHO + 0.05 CH4 + 0.1 GH2 + 0.66 CHAR
    # 3) CELLA -> C6H10O5
    # 4) CELL -> 4.45 H2O + 5.45 CHAR + 0.12 GCOH2 + 0.18 GCH2O + 0.25 GCO + 0.125 GH2 + 0.125 H2
    k1 = 1.5e14 * np.exp(-47_000 / (R * T))
    k2 = 2.5e6 * np.exp(-19_100 / (R * T))
    k3 = 3.3 * T * np.exp(-10_000 / (R * T))
    k4 = 9e7 * np.exp(-31_000 / (R * T))

    # Hemicellulose reactions and rate constants
    # 5) GMSW -> 0.70 HCE1 + 0.30 HCE2
    k5 = 1e10 * np.exp(-31_000 / (R * T))

    # species reaction rate equations where r = dc/dt
    # mass fractions associated with each species are also given
    r_CELL = -k1 * cell
    r_CELLA = k1 * cell - k2 * cella - k3 * cella
    r_CH2OHCHO = k2 * cella * 0.1481
    r_CHOCHO = k2 * cella * 0.0107
    r_CH3CHO = k2 * cella * 0.0462
    r_C6H6O3 = k2 * cella * 0.1944
    r_C2H5CHO = k2 * cella * 0.1254
    r_CH3OH = k2 * cella * 0.0395
    r_CH2O = k2 * cella * 0.0278
    r_CO = k2 * cella * 0.0846
    r_GCO = k2 * cella * 0.008637 + k4 * cell * 0.04319
    r_CO2 = k2 * cella * 0.1167
    r_H2 = k2 * cella * 0.001616 + k4 * cell * 0.001554
    r_H2O = k2 * cella * 0.1033 + k4 * cell * 0.4944
    r_GCH2O = k2 * cella * 0.009259 + k4 * cell * 0.03333
    r_HCOOH = k2 * cella * 0.005677
    r_CH2OHCH2CHO = k2 * cella * 0.02284
    r_CH4 = k2 * cella * 0.004947
    r_GH2 = k2 * cella * 0.001243 + k4 * cell * 0.001554
    r_CHAR = k2 * cella * 0.04889 + k4 * cell * 0.4037
    r_C6H10O5 = k3 * cella
    r_GCOH2 = k4 * cell * 0.02222
    r_GMSW = -k5 * gmsw
    r_HCE1 = k5 * gmsw * 0.7
    r_HCE2 = k5 * gmsw * 0.3

    # system of ODEs
    dcdt = (
        r_CELL, r_CELLA, r_CH2OHCHO, r_CHOCHO, r_CH3CHO, r_C6H6O3, r_C2H5CHO,
        r_CH3OH, r_CH2O, r_CO, r_GCO, r_CO2, r_H2, r_H2O, r_GCH2O, r_HCOOH,
        r_CH2OHCH2CHO, r_CH4, r_GH2, r_CHAR, r_C6H10O5, r_GCOH2, r_GMSW, r_HCE1,
        r_HCE2
    )
    return dcdt

# initial mass fractions [-]
y0 = np.zeros(len(species))
y0[0] = 0.6
y0[22] = 0.4

# solution for system of equations for a batch reactor
sol = solve_ivp(dcdt_debiagi, (time[0], time[-1]), y0, t_eval=time)

# Print
# ----------------------------------------------------------------------------

print('--- Initial ---')
print(f'CELL        {sol.y[0][0]}')
print(f'GMSW        {sol.y[22][0]}')

print('--- Final ---')
for i, sp in enumerate(species):
    print(f'{sp:11} {sol.y[i][-1]:.4f}')
</code>```

Ideally I would like to define all the reactions in an array, parse that array, and output a matrix that can be used by SciPy solve_ivp to determine the concentration of species at each time step. I'm thinking something like what is shown below:

```<code>R = 1.9859      # universal gas constant [cal/(mol K)]
T = 773.15      # temperature [K]

reactions = [
    'CELL -> CELLA',
    'CELLA -> 0.40 CH2OHCHO + 0.03 CHOCHO + 0.17 CH3CHO + 0.25 C6H6O3 + 0.35 C2H5CHO + 0.20 CH3OH + 0.15 CH2O + 0.49 CO + 0.05 GCO + 0.43 CO2 + 0.13 H2 + 0.93 H2O + 0.05 GCH2O + 0.02 HCOOH + 0.05 CH2OHCH2CHO + 0.05 CH4 + 0.1 GH2 + 0.66 CHAR',
    'CELLA -> C6H10O5'
]

names = {
    'CELL': 'C6H10O5',
    'CELLA': 'C6H10O5',
    'GCO': 'CO',
    'GCH2O': 'CH2O',
    'GH2': 'H2',
    'CHAR': 'C'
}

arrhenius = [
    (1.5e14, 0, 47_000),
    (2.5e6, 0, 19_100),
    (3.3, 1, 10_000)
]
</code>```

How can I define a system of chemical equations which can be implemented with solve_ivp without having to write every rate equation by hand for each chemical species?


问:
Normalising DFTs Correctly
I have been playing around with convolutions in scipy's signal package:

```<code>N = 2**16
t = np.linspace(-4, 4, N)
h = rect(t, -1, 1)

scipy_conv = signal.fftconvolve(h, h, mode = 'same')
scipy_conv = scipy_conv/max(scipy_conv)

ax1.plot(t, h, label = r'<span class="math-container">$h(t)$</span>')
ax1.plot(t, scipy_conv, label = r'<span class="math-container">$(h*g)(t)$</span> (SciPy)', linestyle = '--')
ax1.set_xlabel(r'<span class="math-container">$t$</span>')
ax1.grid()
ax1.legend()
</code>```

which as expected produces the plot

However, if I omit the artificial normalising done by diving by the maximum value of the convolution, the convoluted signal has a massive value:

My question is, how should I go about correctly normalising the convolution so that its amplitude would match that of the analytic result of the convolution of the rectangle with itself?

I'm aware this is somehow artefact of the fact we are using DFTs and not continuous FFTs, and that the answer is likely very simple, I am just rather confused. Any help or pointers to an answer I have missed is much appreciated. Thanks in advance.

Title: Normalising DFTs Correctly


问:
What exactly is the cause(s) of blow-up for too-large step size in a method like RK4?
I have been working on creating a few home-made numerical methods, and I am using them to visualize text-book problems from my Strogatz dynamics textbook. It feels like a good way to learn numerical methods and dynamics at the same time.

A recent problem from the book deals with "Critical Slowdown" and is exemplified by comparing $\dot{x} = -x^3$ with $\dot{x} = -x$. I can't pretend to understand anything except that one decays more quickly than the other.

I decided to run both through my little home-made RK4 algorithm and graph them. What I found is that I get overflow for too-large step size, and I was wondering why. Conceptually, what causes this type of behavior?

I know that the step size should be chosen in accordance with how quickly the solution varies, or error will be very high. What I don't get is the overflow errors I'm getting.

Code (The ODE's are defined and solved at the bottom. To see the overflow, replace .001 with .1 in the step-size argument.:

```<code>def rk4(dt, t, field, y_0):

    y = np.asarray(len(t) * [y_0])

    for i in np.arange(len(t) - 1):
        k1 = dt * field(t[i], y[i])
        k2 = dt * field(t[i] + 0.5 * dt, y[i] + 0.5 * k1)
        k3 = dt * field(t[i] + 0.5 * dt, y[i] + 0.5 * k2)
        k4 = dt * field(t[i] + dt, y[i] + k3)
        y[i + 1] = y[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6

    return y

def vf_grapher(fn, t_0, t_n, dt, y_0, lintype='-r', sup_title=None,
               title=None, xlab=None, ylab=None):

    t = np.arange(t_0, t_n, dt)
    y_min = .0
    y_max = .0

    fig, axs = plt.subplots()
    fig.suptitle(sup_title)

    axs.set_title(title)
    axs.set_ylabel(ylab)
    axs.set_xlabel(xlab)

    for iv in np.array(y_0, ndmin=1, copy=False):
        soln = rk4(dt, t, fn, iv)
        plt.plot(t, soln, lintype)
        if y_min > np.min(soln):
            y_min = np.min(soln)
        if y_max < np.max(soln):
            y_max = np.max(soln)

    x = np.linspace(t_0, t_n + dt, 11)
    y = np.linspace(y_min, y_max, 11)

    X, Y = np.meshgrid(x, y)

    theta = np.arctan(fn(X, Y))

    U = np.cos(theta)
    V = np.sin(theta)

    plt.quiver(X, Y, U, V, angles='xy')
    plt.xlim((t_0, t_n - dt))
    plt.ylim((y_min - .05 * y_min, y_max + .05 * y_max))
    plt.show()

if __name__ == '__main__':

    def crit_slow(t, x):
        return -x**3

    def exp_decay(t, y):
        return -y

    # graph an ODE modelling critical slow down (has discontinuities in it's
    second derivative)
    vf_grapher(crit_slow, 0., 10., .001, 10.)

    # graph an ODE modelling exponential decay (compare to critical slowdown)
    vf.vf_grapher(exp_decay, 0., 10., .001, 10.)
</code>```


问:
Direct and Inverse efficient mapping of 3D cartesian positions in a 1D array
I saved a sample of $N$ Cartesian locations $\{x_i, y_j, z_k\}$ inside a one-dimensional array $\mathbf{a} = \{a_l\}_{l = 1}^N$.
How can I access back (efficiently) the $l$-location of the array $\mathbf{a}$ knowing just the three $i,j,k$ Cartesian indices?

(Fortran/C++ examples would be helpful)


问:
Symmetric sparse direct solvers in scipy
scipy.linalg.solve, in its newer versions, has a parameter assume_a that can be used to specify that the matrix $A$ is symmetric or positive definite; in these cases, LDL or Cholesky are used rather than LU (Lapack's sysv and posv rather than gesv).

Is there a similar interface for sparse solvers? As far as I understand, scipy.sparse.linalg.spsolve does not support assume_a and always uses LU. What is the recommended way to use a symmetric sparse direct solver in Scipy, then (if there is any at all)?

I have seen that there is also sksparse.cholmod, but it is a separate package with a different interface, and from the documentation it looks like it does not handle indefinite matrices at all.


问:
Speedup of CPU Pipelining by number of steps
When a CPU has $K$ steps the speed up of using pipelining compared to non-pipelining is $K$. But what I want to know is, say I am a CPU designer and want to decide whether I should build $K$ or $N$ steps, how do I know the speedup between having $K$ or $N$ steps?

I figured that the theoretical speed would be the same and the only possible speedup is when steps can be implemented with simpler circuits, thus increasing the clock speed. Am I right with this, assuming that the two designs only differ by number of steps?

Another point is, how do more steps influence the real-time aspect of processing? Since more steps means waiting long, is there a point, when the pipelining may have a bad impact on the system stability? Bechaus essential what Pipelining does is to value throughput more than speed.


问:
Solving a system of PDEs with no-flux boundary conditions (finite difference discretization)
I am interested in solving a system of linear PDEs with the finite difference method and I'm having trouble to solve the no-flux boundary condition correctly.
\begin{align}
\frac{\partial n}{\partial t} &= D_n\nabla^2n - \nabla\cdot\left(\chi n\nabla c\right) - \nabla\cdot (\rho n\nabla f) \\
\frac{\partial f}{\partial t} &= \beta n-\gamma nf \\
\frac{\partial c}{\partial t} &= -\eta nc
\end{align}
The domain is $\Omega=[0,1]\times[0,1]$ and $D_n ,\chi, \rho, \beta, \gamma, \eta $ are constants.
Also there is the following boundary condition stated, where $\xi$ is the outward unit normal vector:
\begin{align} \left[D_n\nabla n - \left(\chi n\nabla c\right) - (\rho n\nabla f)\right]\cdot \xi = 0
\end{align}.
Considering $\xi = (1,0)$ I received:
\begin{align}
D_n\partial_xn-(\chi n\partial_xc-\rho n\partial_xf)=0
\Leftrightarrow D_n\partial_xn = \chi n\partial_xc-\rho n\partial_xf
\end{align}
So now considering on point on the boundary, eg $(1,n)$ and discreting with finite difference we get:
\begin{align}
N_{1,n+1}-N_{1,n-1} = \frac{\chi}{D_n}N_{1,n}(C_{1,n+1}-C_{1,n-1})-\frac{\rho}{D_n}N_{1,n}(F_{1,n+1}-F_{1,n-1}) \\
\Leftrightarrow N_{1,n+1} = N_{1,n-1} + \frac{\chi}{D_n}N_{1,n}(C_{1,n+1}-C_{1,n-1})-\frac{\rho}{D_n}N_{1,n}(F_{1,n+1}-F_{1,n-1})
\end{align}
So now we got an equation for our ghost point which we can use for equation $n$. The discretization is given as:
\begin{align}
N_{l,m}^{q+1} = N_{l,m}^{q}P_0 + N_{l+1,m}^{q}P_1+N_{l-1,m}^{q}P_2+N_{l,m+1}^{q}P_3+N_{l,m-1}^{q}P_4 \\
F_{l,m}^{q+1} = F_{l,m}^{q}(1-k\gamma N_{l,m}^{q})+k\beta N_{l,m}^{q} \\
C_{l,m}^{q+1} =  C_{l,m}^{q}(1-k\eta N_{l,m}^{q})
\end{align}

if we consider the point $(1,n)$ as stated above we'd receive:
\begin{align}
N_{1,n}^{q+1} = N_{1,n}^{q}P_0 + N_{2,n}^{q}P_1+N_{0,n}^{q}P_2+N_{1,n+1}^{q}P_3+N_{1,n-1}^{q}P_4 \end{align}
The index $(1,n+1)$ is not in our domain, but we can consider the equation for $N_{1,n+1}$ from above and insert.
I'm having trouble when considering the $P_i$, for example:
\begin{align}
P_3 = \frac{kD}{h^2}- \frac{k}{4h^2}(\chi(C_{l,m+1}^{q} - C_{l,m-1}^{q})+\rho(F_{l,m+1}^{q}-F_{l,m-1}^{q}))
\end{align}
My problem is that I don't know how to determine $C_{1,n+1}$ and $F_{1,n+1}$ to get an equation for $N_{1,n+1}$ I can use.

Any suggestions what to try? Please let me know if my question lacks some information, it's my first question here.
Any advice is appreciated!


问:
What is appropriate boundary condition for Poisson pressure equation?
I'm doing CFD simulations in unstructured grids. Well, it's a bit different from conventional unstructured grids that are used mainly in FEM or FVM as tetrahedral meshes. Mine is a voxelized mesh of an unstructured surface, which is shown here as an example:

The left image is my unstructured mesh and the right one is the zoom to show the voxels.

As part of a verification for my CFD simulations, which is written based on lattice Boltzmann method (LBM), I want to check the pressure, which I got from ideal gas equation in LBM as:

$$\Delta P = \Delta \rho c_{s}^{2}$$

Where, $\Delta P$ is pressure difference ($\Delta P = P-P_{ref}$), $\Delta \rho$ is density difference ($\Delta \rho = \rho - \rho_{f}$), $\rho_{f}$ is known fluid density, and finally $c_{s}$ is lattice sound wave velocity ($c_{s} = \hat{c}_{s} \frac{\Delta x}{\Delta t}$, $\hat{c}_{s} = \frac{1}{\sqrt{3}}$, $\Delta x$ is mesh size and $\Delta t$ is time step).

Conventional way of calculating pressure in FEM or FVM is based on Poisson pressure equation:

$$\nabla^{2} P = - \rho_{f}(\nabla \mathbf{u})^{T} : \nabla \mathbf{u}$$

Where $\mathbf{u}$ is velocity of the fluid. I already know the velocity and pressure from LBM, but I want to use just the velocity and put it in Poisson pressure equation to calculate the pressure and compare it to ideal gas equation from LBM. One of the natural Neumann boundary conditions of this Poisson pressure equation is formulated as:

$$\nabla P \cdot \mathbf{n} = \mu \nabla^{2} \mathbf{u} \cdot \mathbf{n}$$

I know this boundary condition is derived from no-slip boundary condition ($\mathbf{u} = 0$) at the boundary $\partial \Omega$. Also, I know at the outlets, I have the pressures. For example, in the simplest case, I put zero pressure boundary condition at the all outlets by assuming the fact that reference pressure is already considered zero: $P = 0$ at $\partial \Omega_{i}^{outlet}$, where $\partial \Omega_{i}^{outlet}$ is the boundary of $i$th outlet. The thing, which I don't know is the pressure boundary condition at the inlet. Initially, in my LBM simulations, I put velocity boundary condition at the inlet and then we could calculate pressure at the inlet plane after simulation is completed. But for extracting pressure from Poisson boundary condition, I don't know what to do at the inlet. I don't want to put calculated pressure from LBM at the inlet cause that somehow will create leakage of pressure information from LBM to this Poisson pressure equation scheme. I can summarize my problem as:

Main equation defined in compuational domain ($\Omega$): $$\nabla^{2} P = - \rho_{f}(\nabla \mathbf{u})^{T} : \nabla \mathbf{u}$$

Boundary conditions:

At the wall ($\partial \Omega$): $$\nabla P \cdot \mathbf{n} = \mu \nabla^{2} \mathbf{u} \cdot \mathbf{n}$$

At the outlets ($\partial \Omega_{i}^{outlet}$): $$P = 0$$

At the inlets ($\partial \Omega_{i}^{inlet}$): ?

Update:

I have an idea for finding a pressure boundary condition at the inlet plane. We have Navier-Stokes as:

$$\rho_{f} \frac{\partial \mathbf{u}}{\partial t} + \rho_{f}\mathbf{u} \cdot \nabla \mathbf{u} = -\nabla P + \mu \nabla^{2} \mathbf{u}$$

If we evaluate Navier-Stokes at the inlet boundary and by assuming the fact that at the inlet plane I put Poiseuille flow, we have:

$$\rho_{f} \frac{\partial \mathbf{u}}{\partial t}|_{\partial \Omega_{i}^{inlet}} \cdot \mathbf{n}_{i}^{inlet} + \rho_{f} \mathbf{u} \cdot \nabla \mathbf{u} |_{\partial \Omega_{i}^{inlet}} \cdot \mathbf{n}_{i}^{inlet} = -\nabla P |_{\partial \Omega_{i}^{inlet}} \cdot \mathbf{n}_{i}^{inlet} + \mu \nabla^{2} \mathbf{u} |_{\partial \Omega_{i}^{inlet}} \cdot \mathbf{n}_{i}^{inlet}$$

Also from Poiseuille flow:

$$\mathbf{u}(\mathbf{r},t)|_{\partial \Omega_{i}^{inlet}} = u_{max}^{i}(t) \Bigg (1 - \frac{|\mathbf{r}-\mathbf{r}_{C}^{i}|^{2}}{R_{i}^{2}} \Bigg ) \mathbf{n}_{i}^{inlet}$$

As a result:

$$\frac{\partial \mathbf{u}}{\partial t}|_{\partial \Omega_{i}^{inlet}} \cdot \mathbf{n}_{i}^{inlet} = \frac{d u_{max}^{i}(t)}{dt} \Bigg (1 - \frac{|\mathbf{r}-\mathbf{r}_{C}^{i}|^{2}}{R_{i}^{2}} \Bigg )$$

$$\mathbf{u} \cdot \nabla \mathbf{u} |_{\partial \Omega_{i}^{inlet}} \cdot \mathbf{n}_{i}^{inlet} = 0$$

$$\mu \nabla^{2} \mathbf{u} |_{\partial \Omega_{i}^{inlet}} \cdot \mathbf{n}_{i}^{inlet} = -\frac{4 \mu u_{max}^{i}(t)}{R_{i}^{2}}$$

Finally:

$$\nabla P |_{\partial \Omega_{i}^{inlet}} \cdot \mathbf{n}_{i}^{inlet} = -\frac{4 \mu u_{max}^{i}(t)}{R_{i}^{2}} - \rho_{f} \frac{d u_{max}^{i}(t)}{d t} \Bigg (1 - \frac{|\mathbf{r}-\mathbf{r}_{C}^{i}|^{2}}{R_{i}^{2}} \Bigg )$$

I appreciate if someone could look at it see if it makes sense or not.


问:
Lax-wendroff for stiff source terms
I am interested in problems of the form
$$
u_t = F(u) + S(u)
$$
where $F(u) = - div(f(u))$ and $S(u)$ is a stiff source term. I am looking for any existing works which develop Lax-Wendroff type schemes for such stiff problems using an IMEX approach to deal with stiff source. So far my search has yielded nothing more than Crank-Nicholson for the source [1] which is not good for stiff case. Is there any fundamental difficulty in constructing Lax-Wendroff schemes which are L-stable ?

The LW schemes I am looking at are of the form as used in [2]. Perhaps LW is not a good name, they seem to be called two-derivative schemes. At second order, they are basically an LW scheme.

[1] Yuangao Zhang and Behrouz Tabarrok, Modifications to the Lax–Wendroff scheme for hyperbolic systems with source terms, https://doi.org/10.1002/(SICI)1097-0207(19990110)44:1%3C27::AID-NME485%3E3.0.CO;2-0

[2] Christlieb et al., Explicit Strong Stability Preserving Multistage Two-Derivative Time-Stepping Schemes, DOI: 10.1007/s10915-016-0164-2

Title: Lax-wendroff for stiff source terms


问:
Why does the matlab command **chol(A)** slower than **chol(A,'lower')** for a large sparse SPD matrix?
For a SPD matrix A, there exists Cholesky factorization $A=LL^T$ or $A=R^TR$, where L, R are a lower and upper triangular matrix, respectively.

Also in matlab, there has a command R = chol(A) which produces $A=R^TR$. and  another command L = chol(A,'lower') which produces $A=LL^T$. But when I implement these two commands with a same large sparse SPD matrix A, L = chol(A,'lower') is faster than R = chol(A) so mysterious. why this happens? Thanks.

clc;clear;
n =400;
A = gallery('poisson',n);

tic
R = chol(A);%   generate the triangular matrix such that A = R'*R
toc

tic
L = chol(A,'lower');%   generate the lower matrix such that A = L*L'
toc
```

<p>I have run the example 3 times and
the numerical results are as follows which indeed demonstrates that the above words what I said (my cpu is 8GB memory and matlab 2018b):</p>

<pre class="lang-matlab prettyprint-override"><code>
Elapsed time is 12.089711 seconds.
Elapsed time is 8.467380 seconds.

Elapsed time is 10.372768 seconds.
Elapsed time is 8.131158 seconds.

Elapsed time is 10.027861 seconds.
Elapsed time is 8.105706 seconds.
</code>```


问:
understanding Domain Decomposition with example
I am new in Domain Decomposition method. I am started to solve $-\Delta u = f$ in $\Omega$ and $u = 0$ on $\partial\Omega$.

From that I get in $\Omega _1$

$$\begin{bmatrix}4&-1\\-1&4\end{bmatrix} \begin{bmatrix}u_1 ^{k +\frac{1}{2}}\\u_2 ^{k +\frac{1}{2}}\end{bmatrix} = \begin{bmatrix}f_1-u_3 ^k\\f_2-u_4 ^k\end{bmatrix}$$

In $\Omega 2$
$$\begin{bmatrix}4&-1\\-1&4\end{bmatrix} \begin{bmatrix}u_3 ^{k + 1}\\u_4 ^{k + 1}\end{bmatrix} = \begin{bmatrix}f_3-u_1 ^{k +\frac{1}{2}}\\f_4-u_2 ^{k +\frac{1}{2}}\end{bmatrix}$$

But in book consider matrix representation as $A_i = (A_{\Omega_i}, A_{\partial\Omega_i-\Gamma_i}, A_{\Gamma_i})$ and for $A_1$ it is given as

$$\begin{bmatrix}4&-1&0&-1&0&0&-1&0&0&-1&0&0\\-1&4&0&0&-1&0&0&-1&0&0&-1&0\end{bmatrix}$$

I did not understand this. Can anyone please explain explicitly?


问:
Can the standard multigrid performance be used for time-dependent PDEs?
Consider a time dependent pde(i.e u(x,t)).I know when only space-coarsening is used the standard multigrid performance can be applied but what if instead we use only time-coarsening?Can we apply the standard multigrid for solving the spatial problem in order to take one time step(i.e going from t1 to t2)?</p>


问:
FEniCS, refinement not 'respecting' domain boundary
Short question: how to ensure that extra points are not included as 'boundary' points after calling the refine function.

More details.

I am working with a hexahedral mesh in $3$d. Let $X$ be the set of mesh points. The mesh describes a connected object with an irregular geometry, but is flat in the sense that if $\pmb{x}\in X$, then the third component of $\pmb{x}$ is in the set $\lbrace z_{1},z_{2}\rbrace$.

The hexahedral mesh was unsuitable for use in FEniCS, so I used the tetrahedralize filter in paraview. I then converted to xml format and was able to read this into a FEniCS code.

The mesh was unsuitable to solve PDEs on since all mesh points were boundary points. I then called refine and have mesh points $\pmb{x}$ such that their z-component is in the set $\lbrace z_{1},z_{3/2},z_{2}\rbrace$. I extracted the boundary mesh by calling bmesh = BoundaryMesh(mesh, "exterior", True). I find that there are too many bmesh coordinate points and that they lie on the interior of the object that I am meshing. Unfortunately the geometry of the mesh is complicated and so I can't describe the boundary with a simple function like $x[0]=0$.


问:
Why wall shear stress calculated from LBM directly and the one calculated based on velocity profile are so different in some cases?
First of all, I hope you accept my apologizes if my question seems off topic here. But, I asked this question in ParaView forum and after a week still I did not receive any response yet, so I'm wondering if I could find an answer hopefully here.

Background: I'm using lattice Boltzmann method to do simulation of blood flow in brain arteries. Basically, those arteries are unstructured manifolds, but we approximate them by using some tiny voxels in order to have an unstructured regular grid (I know it's a bit contradictory to have an unstructured regular grid but please bear with me). I added an example of such voxelized morphology here:

We calculate wall shear stress by using this formula from Matyka et. al., specially equation 2:

$$\vec{\tau}_{WSS} = -\mu \omega \sum_{\alpha} f_{\alpha}^{neq} \vec{c}_{\alpha} \cdot \vec{n} (\vec{c}_{\alpha} - (\vec{c}_{\alpha} \cdot \vec{n}) \vec{n})$$

Where $\vec{\tau}_{WSS}$ is wall shear stress vector, $\mu$ is dynamic viscosity, $\tau$ is relaxation time of LBM framework, $f_{\alpha}^{neq} = f_{\alpha} - f_{\alpha}^{eq}$ is non-equilibrium part of particle distribution, $\vec{c}_{\alpha}$ is discrete velocity in $\alpha$th direction, and $\vec{n}$ is normal vector.

The way initially during voxelization we calculate normal vector of this staircase surface is that: During voxelization, we encounter the intersection of true boundary by using CGAL intersection method, then we assign the intersected normal vector of true boundary to voxels that have intersection with that cell (triangular cell of a STL file if you want to think of it). Then, based on our distance to true boundary, we take an average in each voxel of voxelized surface based on its neighbors and their normal vectors to obtain final normal vector this staircase surface:

$$\vec{n}^{i} = \frac{\sum_{\nu = 1}^{N_{neighbors}} \frac{\vec{n}_{\nu}^{i}}{d_{\nu}}}{\sum_{\nu=1}^{N_{neighbors}}\frac{1}{d_{\nu}}}$$

Where $N_{neighbors}$ are the number of neighbors of the voxel in the surface, $\vec{n}_{\nu}^{i}$ are the initial normal vectors of voxel $i$th $\nu$th neighbor from CGAL, and finally $d_{\nu}$ is the distance of voxel to the true boundary in non-dimensional unit.

So, we get the wall shear stress from here and it is compared to other finite element flow solvers, which shows that error is less than 10%.

On the other hand, from macroscopic point of view, we could calculate wall shear stress from velocity profile as:

$$\vec{\tau}_{WSS} = \mathbf{T} \cdot \vec{n} - (\vec{n} \cdot \mathbf{T} \cdot \vec{n}) \vec{n}$$

Where:

$$\mathbf{T} = \mu (\nabla \mathbf{u} + (\nabla \mathbf{u})^{T})$$

We expect that if we use above formula, we should get more or less same value as LBM calculations. Note that we already know $\mathbf{u}$ from LBM simulation, but we want to use GradientOfUnstructuredDataSet class from ParaView to calculate wall shear stress from velocity field. We do this for 69 different brain artery geometries. Note that we don't use exact same $\vec{n}$ value from voxelization part here. Instead of that, we extract boundary of voxelized mesh by taking ExtractSurface class of ParaView and we use GenerateSurfaceNormals class to create normal vector of that staircase surface.

Problem:

You can see the result of these 69 different cases in this plot, where x axis is the WSS calculated from LBM and y axis is the WSS calculated from ParaView. The thing, which is really strange is that if the above procedure were completely wrong, we should see a systematic offset of datapoints, but obviously my datapoints are all over the place and seems distributed somewhat randomly:

I took, two extreme cases where I have perfect match between LBM and ParaView and where I have complete disagreement between LBM and ParaView and their WSS distribution are shown as below.

Complete disagreement cases:

LBM:

ParaView:

Complete agreement cases:

LBM:

ParaView:

Question:

I appreciate if someone has any idea what's going here. I don't have any idea why it behaves like this. For some cases, it just shows perfect agreement, but for some other cases with exact same script and code it just produce a result completely different. Am I overlooking something hidden in the GradientOfUnstructuredDataSet? Particularly, does anyone know how GradientOfUnstructuredDataSet handles boundaries to take gradient? Still, the answer of those question does not give the answer why for some cases it works quite well and for some other ones it does a pretty bad job.


问:
"Geometry of ill-conditioning" for least-squares problems
It is an idea that dates back to Demmel, 1987 that the condition number of a problem is often related to the distance to the closest ill-posed problems. In Section 3 of the above paper, the author cites as examples linear systems/matrix inversion (1/condition number = relative distance to the closest singular matrix), polynomial root-finding, and eigenvalue calculation (1/condition number = relative distance to the closest problem with double roots).

Can the condition number of least-squares problems $\min \|Ax-b\|$ be interpreted in the same sense?

The first idea that comes to mind is considering 'ill-posed problems' as those with solution $x=0$ and those in which $A$ does not have full column rank, but in this case I cannot make much sense of the $\kappa(A)^2$ that appears in the formulas.


问:
Multilevel minimization - boundary conditions
I am interested in minimizing
$$min_{x \in R^{n^l}} f^l(x),$$
where $f^l(x)$ is nonlinear objective function arising from discretization of PDE. I would like to use nonlinear multilevel minimization techniques, such as MG-opt (Nash, A multigrid approach to discretized optimization problems, https://www.tandfonline.com/doi/full/10.1080/10556780008805795?scroll=top&needAccess=true ).

MG-opt is generalization of full approximation scheme to minimization problems. Hence, on coarse level, it minimizes
$$min_{x \in R^{n^L}} h^L(x) = f^L(x) + \langle R \nabla f^l(x^l_k) - \nabla f^L(R x^l_k) \rangle$$
in order to obtain coarse grid correction.
This objective function is built as a sum of rediscretized fine level physics ($f^L$) and first order consistency term (difference between restricted fine level gradient and initial coarse level gradient).

What is correct way to treat boundary conditions, such that minimization on coarse level is feasible (respects coarse level BC), but does not hinders the convergence properties of the method.  Should coarse level objective function be computed over whole domain, including boundary conditions of the original physics?


问:
Simulating flow in a branched pipe
I am trying to simulate 1D advection and convection of a solute in the following blood vessel segment.
I would like to know if this system can be simulated in COMSOL or MATLAB.

I have used pdepe solver in MATLAB for simulating flows in 1D pipes without branches. However, I am not aware of a tool that can be used for simulating simple branched segments like the above.

Any suggestions?

Title: Simulating flow in a branched pipe