问:
B-spline abscissa values for numerical stability
I am working with spline functions using the B-spline basis,
$$
f(x) = \sum_{i=1}^n c_i B_i(x;\mathbf{t})
$$
where the $B_i$ are cubic B-splines with a knot vector $\mathbf{t}$. For my application, the vector $x$ has units of length, and I am trying to find the best way to nondimensionalize the problem to improve numerical stability. The wiki page on Nondimensionalization suggests simply to scale the values by some appropriate constant,
$$
\chi = \frac{x}{x_s}
$$
where $\chi$ has no units, and $x_s$ is an appropriate scale length for the problem (about $10^6$ for my problem). With this approach, my spline would become,
$$
f(x) = \sum_{i=1}^n c_i B_i(\chi; \mathbf{t}_{\chi})
$$
Using B-splines defined on the $\chi$ variable definitely improves numerical stability when integrating and differentiating the spline. However I am still running into problems. In particular, when computing the 3rd derivative of the spline, I get values on the order of $10^9$ which are causing numerical accuracy problems in later computations.

The matlab algorithm for polynomial fitting, called polyfit, does the following scaling transformation prior to fitting the polynomial,
$$
\chi = \frac{x - \mu}{\sigma}
$$
where $\mu,\sigma$ are the mean and standard deviation of the $x$ vector. The page says that this improves the numerical properties of the fit. I tried using this approach with my spline, and indeed I find that all the derivative terms are close to 1, which seems to be a big improvement.

The purpose of my question is to ask whether the matlab approach is the best one for use with B-splines? The matlab page does not provide a reference to justify its scaling choice, and I have been unsuccessful in finding papers on this topic, especially related to B-splines. Does anyone know of the best way to scale the B-spline knots to get the best numerical accuracy when doing integration/differentation? Does anyone know of a good reference which discusses this problem?


问:
Non-Temporal Weighted Graph Datasets
I am searching for datasets to evaluate an algorithm designed for tasks such as node-classification (edge-prediction, etc.) on weighted and potentially directed graphs.

The  Stanford Network Analysis Project (SNAP) repository only contains datasets containing temporal weighted graphs.

Are there standard datasets for such tasks on non-temporal weighted (and directed) graphs?

ContentLicense: CC BY-SA 4.0

Title: Non-Temporal Weighted Graph Datasets


问:
Navier Stokes Equation Discretization by upwind scheme
I've been reading this paper regarding the SIMPLER algorithm in CFD and I don't understand exactly how the discretized equation mentioned by the author is arrived at. It says that the convective terms are discretized using the upwind scheme and the diffusive terms using central differencing but the final form really doesn't make sense. Also, Reynold's number (Re) is missing in the final form. Below are a few important sections of the paper.

Could someone please walk me through the missing steps?
Thanks in advance.


问:
How do the current FEM opensource libraries compare?
Almost all FEM libraries are good enough, but I want to start with a FEM package and stick to it for some time. Instead of trying all of them, or going with what everyone else is using, I want to benefit from other user's experiences and know:

What are the pros and cons of the FEM libraries you have experience with? How do they compare to other libraries? What are they predominantly used for?

What is the philosophy and vision the developers have?

I already found this link which compares some of the FEM software packages, but I am looking for more input from others.

Other things you can think of when answering this question is, does the FEM package focus on and values:

Generality: the ability to solve almost any problem by providing all the low level tools necessary such as special elements, geometries, special preconditioning, linear solvers...

Efficiency and performance: Focusing on HPC and always trying to squeeze out every bit of computing power one can.

Modernity?: Focusing on the use of state of the art mathematical techniques and hardware.

Low level vs high level: Imagine having two libraries, where it might be easier to use library1 to set up a simple problem, but it can be less efficient than library2 when setting up more complicated problems. This is akin to generality i defined above.

Sorry for keeping things too general.


问:
How to specify the value of pressure in a point in solving steady stokes equation numerically?
I try to solve the steady-Stokes equation numerically on , that is
\begin{aligned}
-\nabla \cdot \mathbb{T}(\mathbf{u}, p) & =\mathbf{f} \quad \text { on } \Omega=[0,1]*[-0.25,0], \\
\nabla \cdot \mathbf{u} & =0 \quad \text { in } \Omega, \\
u_1 & =e^{-y} \text { on } x=0, \\
u_1 & =y^2+e^{-y} \text { on } x=1, \\
u_1 & =\frac{1}{16} x^2+e^{0.25} \text { on } y=-0.25, \\
u_1 & =1 \text { on } y=0 \\
u_2 & =2 \text { on } x=0 \\
u_2 & =-\frac{2}{3} y^3+2 \text { on } x=1, \\
u_2 & =\frac{1}{96} x+2-\pi \sin (\pi x) \text { on } y=-0.25, \\
u_2 & =2-\pi \sin (\pi x) \text { on } y=0 .
\end{aligned}
here \begin{aligned}
& f_1=-2 \nu x^2-2 \nu y^2-\nu e^{-y}+\pi^2 \cos (\pi x) \cos (2 \pi y) \\
& f_2=4 \nu x y-\nu \pi^3 \sin (\pi x)+2 \pi(2-\pi \sin (\pi x)) \sin (2 \pi y)
\end{aligned}
here the $T$ is $$
\mathbb{T}(\mathbf{u}, p)=2 \nu \mathbb{D}(\mathbf{u})-p \mathbb{I}
$$
where $\nu$ is the viscosity and the deformation tensor
$$
\mathbb{D}(\mathbf{u})=\frac{1}{2}\left(\nabla \mathbf{u}+(\nabla \mathbf{u})^t\right)
$$
and the analytical solution is
\begin{aligned}
& u_1=x^2 y^2+e^{-y}, \quad u_2=-\frac{2}{3} x y^3+2-\pi \sin (\pi x) \\
& p=-(2-\pi \sin (\pi x)) \cos (2 \pi y)
\end{aligned}
the procedure of solving this equation numerically is provided by my teacher by using Taylor-hood finite element method,that is

Generate the mesh information matrices $P$ and $T$.

Assemble the stiffness matrix $A$ by using Algorithm I. (We will choose Algorithm I-3 in class)

Assemble the load vector $\vec{b}$ by using Algorithm II. (We will choose Algorithm II-3 in class)

Deal with the Dirichlet boundary condition by using Algorithm III-3.

Fix the pressure at one point in the domain $\Omega$.

Solve $A \vec{X}=\vec{b}$ for $\vec{X}$ by using a direct or iterative method.
after dealing with the Dirichlet boundary condition, I wonder how to fix the value of pressure at a point? I try to fix the value at the first node by eatablishin in the following way:

$A(2*N_b+1,:)=0;\\
A(2*N_b+1,2*N_b+1)=1;\\
b(2*N_b+1,1)=0;\\\\$
here $N_b$ is the number of total basis of quadratic element, is this right?


问:
accelerating solutions of ODEs with close by parameters
Suppose $\mathbf{u}^1\in\mathbb{R}^n$ is an unknown and $\mu^1\in\mathbb{R}$ is a known parameter. Suppose we have solved the non-linear system of equation for $\mathbf{u}^1$ using Newton's iterations.

$$
\mathcal{N}(\mathbf{u}^1;\mu^1)=0
$$

Using $\mathbf{u}^1$ as an initial guess for the Newton iterations, it is possible to accelerate the solution of the system

$$
\mathcal{N}(\mathbf{u}^2;\mu^2) = 0
$$

where $\mu^1$ and $\mu^2$ are sufficiently close. Now, is it possible to extend this to ODEs? For instance, if I have solved

$$
\frac{d\mathbf{u}^1}{dt} = \mathcal{M}(\mathbf{u}^1;\mu^1)
$$

Is it possible to accelerate the solution of

$$
\frac{d\mathbf{u}^2}{dt} = \mathcal{M}(\mathbf{u}^2;\mu^2)
$$

where $\mu^1$ and $\mu^2$ again are sufficiently close.

The only thing I can think of is to integrate both ODEs implicitly with the same time-step. Solve the first ODE, and use the solution at each time-step to accelerate Newton's iterations for the second. Can anything else be done?


问:
Book recommendation on multiphysics
I want to learn multiphysics such as fluid-structure interaction where the simulation is performed for heat transfer, fluid dynamics, and electrodynamics.

Could you please recommend some books about multiphysics that tell you the principles and algorithms?

Title: Book recommendation on multiphysics


问:
Best approach to solve this system of equations?
I have the following 1D (in space, that is) system of equations I would like to solve:

\begin{equation}
    \rho_{fs}\frac{\partial x_{fs}}{\partial t} = h_m\left(W_a - W_{fs}\right) - D_{eff}\left(\frac{\partial \rho_v(T_{fr})}{\partial x} \right)
\end{equation}
\begin{equation}
    \rho_{fr}c_{p,fr}\frac{\partial T}{\partial t} = \frac{\partial}{\partial x} \left( k_{fr} \frac{\partial T}{\partial x} \right) + L_{sub} \frac{\partial \rho_{fr}}{\partial x}
\end{equation}
\begin{equation}
    \frac{\partial \rho_{fr}}{\partial t} = \frac{\partial}{\partial x}\left( D_{eff} \frac{\partial \rho_v(T_{fr})}{\partial x} \right)
\end{equation}

This set of equations describes the time evolution of a frost layer's height, $x_{fs}$, (local) temperature, $T_{fr}$, and (local) density, $\rho_{fr}$.

I have been struggling with trying to solve this system. I have constructed a system of ODEs following the method of lines and attempted to solve it in Matlab using ode45, but this resulted in matrices filled with NaN values.

Following this, I am now trying to solve it through a 1st order upwind scheme for the time derivatives, solving one equation at a time and trying to iterate until the differences are small enough:

\begin{equation}
    \frac{\left(x_{fs}^{n+1} - x_{fs}^n\right)}{\Delta t} = \frac{1}{\rho_{fs}}\left(h_m\left(W_a - W_{fs}\right) - D_{eff}\frac{\rho_v\left(T^n_{i+1}\right) - 2\rho_v\left(T^n_{i}\right) + \rho_v\left(T^n_{i-1}\right)}{(\Delta x)^2} \right)
\end{equation}
\begin{equation}
    \frac{\left(T_i^{n+1} - T_i^n\right)}{\Delta t} = \frac{1}{\rho_{fr,i}c_{p,fr}}\left( k_{fr}\frac{T^n_{i+1} - 2T^n_{i} + T^n_{i-1}}{(\Delta x)^2} + L_{sub}D_{eff}\frac{\rho_v\left(T^n_{i+1}\right) - 2\rho_v\left(T^n_{i}\right) + \rho_v\left(T^n_{i-1}\right)}{(\Delta x)^2} \right)
\end{equation}
\begin{equation}
    \frac{\left(\rho_i^{n+1} - \rho_i^n\right)}{\Delta t} = D_{eff}\frac{\rho_v\left(T^n_{i+1}\right) - 2\rho_v\left(T^n_{i}\right) + \rho_v\left(T^n_{i-1}\right)}{(\Delta x)^2}
\end{equation}

In order to adhere to boundary conditions, the temperature is calculated on the $n_x$ number of nodes, and the density is calculated inbetween, in the $n_x - 1$ number of cells. As for initial and boundary conditions, the temperature at the frost surface and at the bottom of the layer are known, and an initial uniform density profile is assumed due to the very small initial frost height mentioned before.

This approach has not yielded good results yet either, although I am now starting to think it might have to do with my choice of $\Delta t$. The spatial step size $\Delta x$ is set by the layer height $x_{fs}$ and the number of nodes. The starting layer height is very small, in the order of $10^{-4}$. If I understand correctly, this means my $\Delta t$ has to be very small as well in order to satisfy the Courant condition? This would be a bit impractical though, as the simulated time is in the order of tens of minutes or even hours.

I am wondering if there is a better approach to take to solve this system of equations, or if I am making a mistake somewhere. I would appreciate any insights. Apologies if I have made any mistakes in posing my question here.


问:
Finite-difference produces a derivative off by one order
I have a nonlinear 2nd order boundary value differential equation where I used finite-difference method (central finite-difference) to solve it.

$$z''(x)-\frac{\frac{1}{100} z(x)^4 \left(2 z'(x)^2+12\right)-600 \left(z'(x)^2+1\right)-\frac{3 z(x)^8}{500000}}{20 z(x) \left(10-\frac{z(x)^4}{1000}\right)}=0$$

Boundary conditions are,

$$z(0) = 93/10, \quad z(1) = 10^{-3}, \qquad x \in [0,1]$$

My interest is calculating both $z(x)$ and $z'(x)$, I have used $n+1$ points and x[0] as starting value for this calculation. In the end, after obtaining the data points for $z(x)$ I used interpolation to construct the function zfunc to represent $z(x)$. Taking the derivative of zfunc gives zfunc'.

I know that the derivative of $z(x)$ at one of the boundary is $z'(0) \approx -15.53574063$, but as you see zfunc'[0] = -1.50198076. I don't mind if the decimals are a little off but you can see that they are off by an order $O(10^1)$! This is unacceptable, is there something wrong with the finite-difference method that I wrote? I'm guessing it might have something to do with my starting values. I chose those starting values because I already have a rough idea of the curve and my choice is close to that curve.

The plot of zfunc[x] is shown below.

```<code>ClearAll["Global`*"]
Needs["VariationalMethods`"]
f = 1 - (z[x]/zh)^(d + 1);
L = Sqrt[1 + (z'[x]^2/f)]/z[x]^d;(*Lagrangian*)
eulageq = EulerEquations[L, z[x], x];(*Euler-Lagrange equation*)
s = Solve[eulageq, z''[x]][[1]] // Simplify;(*2nd order EOM*)
eq0 = z''[x] - s[[1, 2]] /. {d -> 3, zh -> 10};

a = 0;
b = 1;
n = 250;
h = (b - a)/(n + 1);
alpha = 93/10;
beta = 10^-3;
Z = Array[z, {n + 2}];
z[1] = alpha;
z[n + 2] = beta;
rul = Table[{z''[x] -> ((z[i + 1] - 2 z[i] + z[i - 1])/h^2), z'[x] -> ((z[i + 1] - z[i - 1])/(2 h)), z[x] -> z[i]}, {i, 2, n + 1}];
eqs = Table[eq0 /. rul[[i]], {i, Length[rul]}];
DFx = ResourceFunction["JacobianMatrix"][eqs, Table[z[i], {i, 2, n + 1}]];
x[0] = Table[(1 - i) alpha, {i, 1/100, 90/100, (90/100 - 1/100)/(n)}];

sol = FindRoot[eqs, Table[{z[i], x[0][[i]]}, {i, 2, n + 1}], Jacobian -> DFx, MaxIterations -> 1000, WorkingPrecision -> 30];

xx = Range[a, b, h];
lst = Table[{xx[[i]], z[i] /. sol}, {i, n + 2}];
zfunc = Interpolation[lst, InterpolationOrder -> 3]

zfunc'[0]
-1.5019807629446601256947800

</code>```

I used Mathematica for the code.


问:
Solving PDE on a non-uniform grid with Crank-Nicolson scheme
I am solving a 1D diffusion-type equation with the finite-difference Crank-Nicolson (CN) scheme, and I need to densify the spatial grid around the central point. One could change the spatial variable and re-derive the whole equation in terms of the new variable. Alternatively, in principle it should be possible to do something simpler, that is to re-write all the CN coefficients without assuming a uniform spacing. For example, for the advection term:
$$
    \frac{\partial f}{\partial t} = - \, u\, \frac{\partial f}{\partial x}
$$
this would correspond to replacing
$$
f^{n+1}_i \left[ 1 + u_i \, \frac{\Delta t}{2 \Delta x} \right] - f^{n+1}_{i-1} \cdot u_i \, \frac{\Delta t}{2 \Delta x} = f^{n}_i \left[ 1 - u_i \, \frac{\Delta t}{2 \Delta x} \right] + f^{n}_{i-1} \cdot u_i \, \frac{\Delta t}{2 \Delta x}
$$
with the expression:
$$
f^{n+1}_i \left[ 1 + u_i \, \frac{\Delta t}{2 (x_i - x_{i-1})} \right] - f^{n+1}_{i-1} \cdot u_i \, \frac{\Delta t}{2 (x_i - x_{i-1})} = f^{n}_i \left[ 1 - u_i \, \frac{\Delta t}{2 (x_i - x_{i-1})} \right] + f^{n}_{i-1} \cdot u_i \, \frac{\Delta t}{2 (x_i - x_{i-1})},
$$
where $n$ is the time index and $i$ is the space index. The diffusion term changes accordingly, although with a longer and a bit more complicated expression.

This way, it should be easy to control the spacing as needed based on physical reasons, according to functions that are as smooth as possible. However, this approach does not seem to work. Is there any reason for that? Or some specific constraints that are necessary?

Update: considering the grid built with the spacing in figure 1

I obtain the following solutions, at different energies (figure 2), from where it shows that the low-energy solutions have problems:


问:
Solving PDE on a non-uniform grid with Crank-Nicolson scheme
I am solving a 1D diffusion-type equation with the finite-difference Crank-Nicolson (CN) scheme, and I need to densify the spatial grid around the central point. One could change the spatial variable and re-derive the whole equation in terms of the new variable. Alternatively, in principle it should be possible to do something simpler, that is to re-write all the CN coefficients without assuming a uniform spacing. For example, for the advection term:
$$
    \frac{\partial f}{\partial t} = - \, u\, \frac{\partial f}{\partial x}
$$
this would correspond to replacing
$$
f^{n+1}_i \left[ 1 + u_i \, \frac{\Delta t}{2 \Delta x} \right] - f^{n+1}_{i-1} \cdot u_i \, \frac{\Delta t}{2 \Delta x} = f^{n}_i \left[ 1 - u_i \, \frac{\Delta t}{2 \Delta x} \right] + f^{n}_{i-1} \cdot u_i \, \frac{\Delta t}{2 \Delta x}
$$
with the expression:
$$
f^{n+1}_i \left[ 1 + u_i \, \frac{\Delta t}{2 (x_i - x_{i-1})} \right] - f^{n+1}_{i-1} \cdot u_i \, \frac{\Delta t}{2 (x_i - x_{i-1})} = f^{n}_i \left[ 1 - u_i \, \frac{\Delta t}{2 (x_i - x_{i-1})} \right] + f^{n}_{i-1} \cdot u_i \, \frac{\Delta t}{2 (x_i - x_{i-1})},
$$
where $n$ is the time index and $i$ is the space index. The diffusion term changes accordingly, although with a longer and a bit more complicated expression.

This way, it should be easy to control the spacing as needed based on physical reasons, according to functions that are as smooth as possible. However, this approach does not seem to work. Is there any reason for that? Or some specific constraints that are necessary?

Thanks a lot!

OwnerDisplayName: ottavio


问:
FEM Basis for functions of two variables in $\mathbb{R}^2$ (Applied to Linear Full Stokes Equations)
I'm trying to do FEM for a very basic version of the linear full Stokes equations in two dimensions. Say we are working in the grid $[0,1]\times[0,1]$ in the $xy-$plane.

To solve an FEM problem for a scalar function $u(x,y)$ in the past, I have used the 2D FEM basis function definition $\psi_{ab}(x,y) = \phi_{a}(x)\phi_{b}(y)$ for $1\leq a,b \leq N-1$ with

$$
\begin{equation*}
    {\phi _a (x) =\left\{
        \begin{array}{ll}
          \frac{x-x_{a-1}}{h} & x\in [(a-1)h,ah] \\
          \frac{x_{a+1}-x}{h} & x \in [ah,(a+1)h] \\
          0 & \mbox{else}
        \end{array}
    \right.}
    \mbox{   }\mbox{   } \mbox{   }\mbox{   }\mbox{   }\mbox{   for } a = 1,...,N-1
\end{equation*}
$$

$N = $ number of intervals we cut the full interval $[0,1]$ into (can change for $x$ and $y$ intervals)

$h = 1/N$

$x_{j} = jh$

I now want to create a basis to represent a vector function $v(x,y)\in \mathbb{R}^2$ for the 2D linear full stokes equations. My current thinking is that if I have some vector function

$$
\begin{equation*}
    v (x,y) =\left(
        \begin{array}{ll}
          v^1 (x,y) \\
          v^2(x,y)
        \end{array}
    \right)
\end{equation*}
$$

I think I can represent this function as $ v = \sum _ {a,b = 1} ^ {N-1} V_{a,b} \hat{\psi}_{ab}(x,y)$ where $V_{a,b}$ are coefficients and

$$
\begin{equation*}
    \hat{\psi}_{ab}(x,y) =\left(
        \begin{array}{ll}
          \phi_{a}^1(x)\phi_{b}^1(y) \\
          \phi_{a}^2(x)\phi_{b}^2(y)
        \end{array}
    \right)
\end{equation*}
$$

My thinking is if we know that $\psi(x,y)$ is a basis (since I was given the first definition I wrote), then we know that
$c_{11} \psi_{11} + c_{12} \psi_{12}+...+c_{N-1,N-1} \psi_{N-1,N-1} = 0 $

if and only if

$$ c_{11} = c_{12} = c_{13}=...=c_{N-1,N-1}=0$$

So I would think that if you took a linear combination of the new basis I wrote for the vector-valued basis functions we would get
$$c_{11} \hat{\psi}_{11} + c_{12} \hat{\psi}_{12}+...+c_{N-1,N-1} \hat{\psi}_{N-1,N-1} = \vec{0} $$ if and only if

$$ c_{11} = c_{12} = c_{13}=...=c_{N-1,N-1}= \left(\begin{array}{ll}
          0 \\
          0
        \end{array}\right) $$

Which I think shows linear independence.

Is what I outlined the correct FEM basis to represent a function in
$\mathbb{R}^2$?

How would I even interpret the coefficients $V_{ab}$ here? Previously, in the case when I used the $\psi(x,y)$ basis functions $V_{ab}$ was a matrix. Now we get a coefficient that is a 2D vector for each combination of a and b (instead of a single number). How can I visualize this? Do I just input a vector into a matrix each time to track the coefficients?

Sorry if these are silly questions I'm struggling to fully visualize this and can't seem to find anything very helpful online.


Title: Solving the TOV equation using Python code: How to manage with the pressure at the surface?

Body:
I wrote the code to solve the TOV equation using the Euler method but I got a negative density and pressure! I could not fix the while loop correctly (steps until reaching the surface.) Can anyone tell me how to fix this error? This is my code:</p>
```<code>#============================================================================
def Euler(t, ht, y, n, Func):
#----------------------------------------------------------------------------
#  Propagates the solution y of a system of 1st order ODEs
#     y'[i] = f[i](t,y[]), i = 1..n
#  from t to t+ht using Euler's method
#  Calls: Func(t, y, f) - RHS of ODEs
#----------------------------------------------------------------------------
   f = [0]*(n+1)                                                # RHS of ODEs

   Func(t,y,f)                                              # get RHS of ODEs
   for i in range(1,n+1): y[i] += ht * f[i]              # propagate solution
   return t, y
#============================================================================

#----------------------------------------------------------------------------
#  Solving the TOV equations by Euler's method
#----------------------------------------------------------------------------

from math import pi

# Some constants
G = 6.67259 * pow(10., -8.)           # gravitational constant  (dyne.cm**2/s)
M_sun = 1.98892 * pow(10., 33.)       # mass of the Sun   ( g )
c = 2.99792458 * pow(10., 10.)        # speed of light    (cm/s)
km = 1. * pow(10., 5.)                # km in cm's

# The following are the quantities at the origin r = 0
k = pow(10., 10.)
p_c = 1.7*pow(10., 35.)
rho_c = (p_c/k)**(3./5.)

# Near the origin, these quantities, and others take the form
delta_r = 100.
r0 = delta_r
m0 = (4./3.)*pi*rho_c*r0**3
xi = (1 - (p_c/(rho_c*c**2)))*(1 + (3.*p_c/(rho_c*c**2)))
p0 = p_c - (2.*pi/3.)*G*rho_c**2*r0**2*xi
rho0 = (p0/k)**(3./5.)

print ('r_ini        mass_ini       pressure_ini\n')
print('{} {} {}'.format(r0/km, m0/M_sun, p0))

"""Now the RHS of the ODEs where y(r) refers to p(r) and m(r)
Analytical expressions (from Shapiro, Teukolsky, Eqts (5.7.5), (5.7.6)):
dm/dr = 4.*pi*rho*r**2 equiv to y[1]' = 4.*pi*rho*r**2
dp/dr = - (rho*m/r**2)*(1 + p/rho)*(1 + (4.*pi*p*r**3)/m)*(r/(r - 2.*m))
equiv to
y[2]' = - (rho*y[1]/r**2)*(1 + y[2]/rho)*(1 +\
 (4.*pi*r**3*y[2])/y[1])*(r/(r - 2.*y[1]))"""
# Surface pressure
p_s = pow(10., 6.)
def Func(r, y, f):
    m = y[1]
    p = y[2]
    if p < p_s:
        p = 0.0
        rho = 0.0
    rho = (p/k)**(3./5.)
    f[1] = 4.*pi*rho*r**2
    if p < p_s:
        f[2] = 0.0
    else:
        f[2] = - G*(rho*m/r**2.)*(1 + p/(rho*c**2.))\
        *(1. + (4.*pi*r**3*p)/(m*c**2.))*(r*c**2./(r*c**2. - 2.*G*m))
    return f

dr = 100.                                   # step size
r = dr
n = 2                                      # number of 1st order ODEs
y = [0]*(n+1)                              # solution components

out = open("NstarGR.txt","w")             # open output file
out.write("#r            Mass      density     Pressure\n")

y[1] = m0; y[2] = p0                         # initial values
out.write(('{} {} {} {}\n').format(r/km, y[1]/M_sun, rho_c, y[2]))

# Now, loop until the surface
p = y[2]
while (p > p_s):
    Euler(r, dr, y, n, Func)
    p = y[2]
    rho = (p/k)**(3./5.)
    r += dr
    out.write(('{} {} {} {}\n').format(r/km, y[1]/M_sun, rho, p))
out.close()
print ('R        M       P\n')
print('{} {} {}'.format(r/km,y[1]/M_sun, p))
</code>```


问:
(Algorithmic) Differentiation capable Finite Volume Software: Generation Jacobian
I am looking for Finite Volume Software that employs a method-of-lines like approach by constructing from the hyperbolic PDE of form
$$\partial_t \boldsymbol u(t,\boldsymbol x) + \nabla \cdot \boldsymbol f\big(\boldsymbol u(t,\boldsymbol x)\big) = \boldsymbol 0 \tag{1} $$
the semidiscretized ODE
$$ \frac{d }{d t} \boldsymbol U(t) =  \boldsymbol F\big(\boldsymbol U(t)\big) \tag{2} \label{2} $$
and allows the computation of the Jacobian of the RHS of \eqref{2}:

$$ J(\boldsymbol u) = \nabla_{\boldsymbol u} \boldsymbol F (\boldsymbol u) $$

Preferably, this is done with algorithmic/automatic differentiation, although finite differences should also be alright for the beginning.
Also, computation for the initial value $\boldsymbol u_0$ would also suffice here.

I checked standard software like Clawpack and FiPy but they do not seem to offer such capabilities.


问:
Open source implementations of the medial axis transform for vector shapes
Are there any open source implementations of the medial axis transform for vector shapes?

I have searched without finding any useful results. It seems that CGAL library doesn't have it implemented nor the deal.ii library.

Are there any impelementation?

Note: Some examples of shapes I am interested in: Polygons with holes. The holes are either circular or polygonal.


问:
Sampling points from space based on known density
I have a problem where I heavily need to restrict the number of points at which I sample a function based on the values of a different function.

I have two functions:

$f:{\mathbb{R}\times [0,\infty)\times[0,\infty) }\rightarrow\mathbb{R}$, which is a probability density, generated on a set of data, using Gaussian kernels

$\sigma:\mathbb{R}^3\rightarrow\mathbb{R}$, which contains information about the quality of my model.

I try to evaluate to total quality of my model, by sampling points from the space $[a,b] \times [c,d] \times[e,f]$ and integrating over the product $f \cdot \sigma$. The intervals correspond to the natural boundaries of my data.

$f$ has at least two peaks, but could contain more interesting points, yet $f([a,b] \times [c,d] \times[e,f])$ evaluates to $0$ in most points.
Now to the problem:

$f$ is relativeley easy to compute for a big amount of points, while $\sigma$ is computationally quite heavy. Therefore I try to find a way to sample $\sigma$ (or both functions) at fewer points, by using $f$ as a literal "density", sampling more points in regions of local extremes, or a higher probability density, while sampling fewer points in the regions where $f$ is $0$ or near $0$ anyway.

Is there a good/efficient way to do so? Main goal is avoiding skipping over intersting regions in $f$ while keeping the number of samples as low as possible. Bonus points: Is there MATLAB funtionality to do so, already?


问:
Does the choice of a complex inner product affect Krylov methods?
As far as I understand there are two definitions of the complex inner product:
$$(a,b) = b^H a$$
and
$$(a,b) = a^H b$$
I know some linear algebra libraries such as BLAS and Eigen uses the second one.

However, I've seen some papers, such as this one, that use the first one.
I would like to know if using one instead of the other would dramatically affect the behavior of a Krylov method such as the CG or the BiCGSTAB.


问:
(Dynamic FE analysis) Proportionality between time step and displacement
I am currently adopting a Newmark algorithm, for which the system of equation is

$$(M+{\beta}{\Delta}t^2{K}){d_i}=\beta{\Delta}t^2{F_i}+Md_{i-1}+{\Delta}tM\dot{d}_{i-1}+{\Delta}t^2M\left(\frac{1}{2}-\beta\right)\ddot{d}_{i-1}$$

Say I solve the problem for a given time step ${\Delta}t$ and find $d_{i}$. If I want to find ${\alpha}d_{i}$ where $\alpha$ is known, by which amount must the time step be multiplied?


问:
Parameter choice rules for L1 regularization?
I am solving an L1 regularized least squares of the form like:

$$ \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| A \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda {\left\| \boldsymbol{x} \right\|}_{1} $$

I saw that for the L2 norm case there are several methods to obtain the magnitude of $\lambda$ as seen in the study Comparing parameter choice methods for regularization of ill-posed problems. For the L1 case I was told to look at AIC\BIC but I have a feeling I am not using them correctly.

For example, y is a 1x1000 vector and A is a 1000x1000 dictionary, and I plan it so that x will be a sparse vector of a length 1000 with only 5-15 non zeros values. By trial and error the range of $\lambda$s that retrieve accurate solutions is broad in the range 1e-4-1e-2 (using CVX).
I want to find the optimal value, similar to generalized cross validation that is done in L2 norm case.

I am calculating AIC\BIC (using aicbic function in matlab) as follows: I scan values of $\lambda$ in the range 1e-7 to 1. For each $\lambda$ value I use for the log-likelihood, the log of the output of the optimal value of CVX (cvx_optval), for the number of estimated model parameters I use the number of non-zero (or |x|>1e-6) entries in the solution of x, and for the sample size I simply use the length of y.

I get the minimum of the AIC\BIC estimate for a $\lambda$ value around 1, which is incorrect. This is because the number of non-zero entries that is used for the estimated model parameters is dropping to 1, so the total function gets a minima there.

What am I doing wrong, and should I use AIC\BIC or something else and how?


问:
Adjoint method for imaging with an "analytical" forward model
I am working on tomographic image reconstruction (radar regime) of the dielectric properties of objects. As part of my work, I have programmed a ray tracer. This ray tracer is, in a way, purely "analytical" in that rays are propagated from their current location to the next intersection with a surface of an object or the receiver (if applicable), always obeying Snell's law when changing the medium and otherwise travelling in straight lines. In other words: this ray tracer does not use a FEM or anything of the sorts; both signal attenuation and optical path length are simply "accumulated" for each individual ray. Importantly, the ray tracer can follow multiple signal propagation paths. All rays that hit the receiver are summed up.

I want to use this ray tracer as a forward model for image reconstruction, and I want to use the full signal for both the simulated and measured data. My first idea was to use the conjugate gradient method. It seemed very practical because you want to minimize a penalty term, e.g.:
\begin{align}
f = \sum_x (s_{\text{sim}}(x) - s_{\text{meas}}(x))^2 \text{,}
\end{align}
where $s_{\text{sim}}(x)$ is the simulated (complex) signal at position $x$ and $s_{\text{meas}}(x)$ the corresponding measured signal. (Note that these signals consist of a frequency chirp, i.e. the signal at each position is a vector.) Minimizing the penalty term should mean that the objects in the ray tracer have become very close to the true objects in terms of their dielectric properties. (We also assume that the position and size of the objects we put into the ray tracer are already correctly determined.)

However, I realized that this method necessitates the computation of the gradient of $f$ and therefore of the simulated signal with regards to the parameters $p$, i.e. $\nabla_p s_{\text{sim}}(x)$. Here, the parameters are $\mu$ and $n$, i.e. the attenuation coefficient and refractive index of the material. This quickly gets very computationally expensive even with a ray tracer as soon as a few objects are used, since for each object and each parameter the gradient needs to be calculated individually (and then also all possible combinations).

I found out that the adjoint method can solve this problem; it can indirectly determine the gradient with regards to the parameters by instead taking the gradient with regards to the simulated signal. I have also read that you need to solve the adjoint equation, and you basically "backpropagate" the residual signal (i.e. the difference between simulated and measured signal). However, most of the applications I found focus on a forward model which uses a FEM. Since that is not the case here, I have been struggling with finding out how to proceed.

Does anyone have suggestions on how I can implement the adjoint method here, or maybe can point me towards literature that might be helpful in my case?


问:
Want to know how to solve this high order differential equation system with boundary condition?

The numerical results may be like this:

Body: The numerical results may be like this:


问:
Numerical Method for Multivariate Inversion Formula
For my research, I need to evaluate the density of a random vector $\boldsymbol{X} \in \mathbb{R}^p$ using the multivariate inversion formula. Let the density function of $\boldsymbol{X}$ be $f (\boldsymbol{x})$. If the characteristic function of $\boldsymbol{X}$ is given by
$$
\phi(\boldsymbol{t}) = \int_{-\infty}^{\infty} \dots \int_{-\infty}^{\infty} e^{i \boldsymbol{t}^\top \boldsymbol{x}} f(\boldsymbol{x}) d\boldsymbol{x},
$$
then using the multivariate formula,
$$
f (\boldsymbol{x}) = (2 \pi)^{-p} \int_{-\infty}^{\infty} \dots \int_{-\infty}^{\infty} e^{-i \boldsymbol{t}^\top \boldsymbol{x}} \phi (\boldsymbol{t}) d\boldsymbol{t}.
$$
In the paper of Sheppard (1991), there is a numerical procedure through the Riemann sums mentioned to evaluate the density function, but it is quite complicated and even gets more so for higher dimensions. Plus, given that I don't have very good knowledge of numerical methods, I am really looking for methods that are already implemented or at least not too hard to implement by myself. My working programming language is R.

I am actually aware of two R packages that can approximate the Fourier integrals: CharFunToolR and fourierin. Both of them are limited to lower $p < 3$: CharFunToolR only works for univariate data while fourierin only allows for up to bivariate now. One problem with the fourierin package is that it requires the upper and lower limits be real values. When I tested with some know distributions, for a given set of upper and lower limits, fourierin did not yield correct densities for all values (usually those in the boundary; the densities are very different from the truth). This is not convenient because the estimated densities may not be reliable.

From my search so far, I understand that this is a hard problem because we are integrating over the infinite domains and higher dimensions are numerically challenging. For that reason, I guess at this point, I am okay with dimensions within 5 (10 would be nice but I think I cannot be too greedy). Please help me if you have any ideas. Thank you so much.

Here is a code snippet regarding the multivariate normal distribution. I put here just in case it is useful somehow.

```<code>p <- 6

mu <- c(-0.448, -0.48, -0.661, -0.325, -0.539, -0.297)

Sigma <- matrix(nrow = p, ncol = p, c(
   0.136, -0.013,  0.001, -0.009, -0.005, 0.004,
  -0.013,  0.123,  0.050,  0.002,  0.037, 0.040,
   0.001,  0.050,  0.095, -0.011,  0.056, 0.014,
  -0.009,  0.002, -0.011,  0.099, -0.005, 0.008,
  -0.005,  0.037,  0.056, -0.005,  0.075, 0.003,
   0.004,  0.040,  0.014,  0.008,  0.003, 0.102
))

x <- c(0.053, 0.523, -0.041, -0.589, -0.044, -0.066)

cf <- function(tee, x, mu, Sigma) {
  exp(-1i * crossprod(tee, x)) * exp(1i * crossprod(tee, mu) - 0.5 * t(tee) %*% Sigma %*% tee)
}

mvtnorm::dmvnorm(x, mu, as.matrix(Sigma))
</code>```


问:
Normalisation in tensor networks
I am trying to implement the iTEBD algorithm for the $PXP$ model, i.e, the hamiltonian is

$$H = \sum_iP_{i-1}X_iP_{i+1}.$$

Here $P$ is the projector onto the ground state and $X$ is the usual pauli x matrix.

I aim to time evolve a certain initial state and record certain observables. To do so, I initialized my iMPS (infinite Matrix Product State), and as a check, I took the overlap of this with itself. As I expected, this gave me a result of 1. However, it appears that as I time evolve my state, the normalization strays further from unity, and by the last time step, it is of the order 1e5. I know that we perform a truncation in each time step, but my question is can this introduce such a large error?  Or is there a specific step where we have to normalize the MPS that I am missing?

Title: Normalisation in tensor networks


问:
ValueError: array must not contain infs or NaNs; When using solve_ivp in the scipy library
I am solving an initial value problem using solve_ivp. The problem consists of computing the concentration profile of a set of reactions over time, given the initial concentrations and some of the reaction rate parameters. Furthermore, I have some experimental data of those concentration profiles at different temperatures for which I have to determine those reaction rate constants by optimizing those parameters. Here is my code implementation:

```<code>import numpy as np
import scipy
import matplotlib.pyplot as plt
import pandas as pd

Npoints = 1001
t = np.linspace(0, 1e-6, Npoints)
IC = np.array([0, 0, 0.5, 1.5, 0]) # initial conditions
species = ['Br', 'H', 'Br2', 'H2', 'HBr']
temperature = [1300, 1400, 1500]
exp_data1300 = pd.read_excel("File", sheet_name="1300K", usecols=[1,2,3,4,5])
exp_data1400 = pd.read_excel("File", sheet_name="1400K", usecols=[1,2,3,4,5])
exp_data1500 = pd.read_excel("File", sheet_name="1500K", usecols=[1,2,3,4,5])
exp_time = pd.read_excel("File", sheet_name="1300K", usecols=[0])
exp_data1300 = exp_data1300.to_numpy()
exp_data1400 = exp_data1400.to_numpy()
exp_data1500 = exp_data1500.to_numpy()
exp_time = exp_time.to_numpy()

def chain_HBr(t, conc, T, params):
    # Exctract concentration data and parameters
    cbr, ch, cbr2, ch2, chbr = conc
    Ea1, A1, m1, Ea2, A2, m2 = params
    R = 8.314 # Gas constant
    A3=5.41e9
    m3=2.05
    Ea3=7550
    A4=5.01e9
    m4=1.05
    Ea4=-680
    A5=5.06e8
    m5=0.1
    Ea5=-8400
    # Define reaction rate constants
    k1 = A1*(T/298)**m1*np.exp(-Ea1/R/T)
    k2 = A2*(T/298)**m2*np.exp(-Ea2/R/T)
    k3 = A3*(T/298)**m3*np.exp(-Ea3/R/T)
    k4 = A4*(T/298)**m4*np.exp(-Ea4/R/T)
    k5 = A5*(T/298)**m5*np.exp(-Ea5/R/T)
    # Define reaction rates of each reaction
    r1 = k1*cbr2
    r2 = k2*cbr*ch2
    r3 = k3*ch*cbr2
    r4 = k4*ch*chbr
    r5 = k5*cbr**2
    # Define species reaction rate
    dBrdt = 2*r1-r2+r3+r4-2*r5
    dHdt = r2-r3-r4
    dBr2dt = -r1-r3+r5
    dH2dt = -r2+r4
    dHBrdt = r2+r3-r4
    return np.array([dBrdt, dHdt, dBr2dt, dH2dt, dHBrdt])

def errorfunc(params, t, c_exp1300, c_exp1400, c_exp1500, y0):
    # Solve concentration profile at each temperature
    sol1300 = scipy.integrate.solve_ivp(chain_HBr,[t[0],t[-1]],y0,t_eval=t,args=(1300,params),method='Radau',max_step=0.01)
    sol1400 = scipy.integrate.solve_ivp(chain_HBr,[t[0],t[-1]],y0,t_eval=t,args=(1400,params),method='Radau',max_step=0.01)
    sol1500 = scipy.integrate.solve_ivp(chain_HBr,[t[0],t[-1]],y0,t_eval=t,args=(1500,params),method='Radau',max_step=0.01)
    # Compute Residual sum of squares
    error1300 = np.sum((sol1300.y.transpose()-c_exp1300)**2)
    error1400 = np.sum((sol1400.y.transpose()-c_exp1400)**2)
    error1500 = np.sum((sol1500.y.transpose()-c_exp1500)**2)
    residual = error1300+error1400+error1500
    return residual

guess = np.array([-2e5, 4e12, -0.5, -7e4, 5e10, 0.5])
params_fit = scipy.optimize.least_squares(errorfunc,guess,args=(exp_time.reshape((len(exp_time))),exp_data1300,exp_data1400,exp_data1500,IC), max_nfev=10000)
</code>```

But most solvers simply return: ValueError: array must not contain infs or NaNs because the time scale is of the order of 10^(-7), which gives a massive differential term in each iteration of solve_ivp. I have tried different solvers like LSODA, but that one just keeps computing an answer forever, and I have also tried setting a smaller max_step, but a max_step of 0.00001 still gives the same error, and anything smaller than that will take ages to solve on my computer. So, does anyone have any other tips in how to go about implementing this? Or what modifications could I make to give me a result?

Thanks!


问:
How to perform a monte carlo simulation on a honeycomb lattice with first, second and third nearest neighbor interactions?
I want to perform monte carlo simulation on a honeycomb lattice. Can someone please help me on how to define a matrix for these interactions and how many terms there will be while sampling the configuration? A honeycomb lattice has 3 first nearest-neighbors, 6 second nearest-neighbors and 3 third nearest-neighbors.</p>


问:
Solving $u_{tt} = u_{xx}$ with spectral methods and homogeneous Dirichlet/Neumann BC
Suppose we want to approximate a function $u(x,t)$ with $(x,t)\in[-1,1]\times[0,4]$ that satisfies the wave equation $u_{tt}=u_{xx}$ in the inside and BC $u(-1,t)=u_x(1,t)=0$ with initial value given by $u(x,0)=\exp(-200x^2)$ and $u_t(x,0)=u_x(x,0)$ (so it is a left-moving Gaussian pulse). The solution would be a pulse that bounces at $T=1$ and changes amplitude to negative values and goes to the right until $T=3$ where it bounces without changing the sign and arrives at $T=4$ at the middle. The "theoretical solution" that we were thinking was to double the domain to the right to get $x\in[-1,3]$ and then the solution would be decomposed into cosines, however, we don't care about that because we want to solve it numerically with Chebyshev spectral methods. So we are using Trefethen's Spectral methods for MATLAB [D,x]=cheb(N) to find a discretization in time, and if it were homogeneous Dirichlet BC's then all we would had to do would be to leapfrog with $\exp(-200x^2)$ and $\exp(-200(x-dt)^2)$ all the way until the end, multiplying the interior points of $x$ by $D$ and ignoring the contribution by the boundary points as we can take them from $0$. We have seen the trick in program 33 of Trefethen when solving a BVP $u''=f$ that they take a row D^2(end,:)=D(end,:) to then solve x=D\f. But for the leapfrog we don't even know where to start, do we have to solve a system at each timestep? How do we enforce the Neumann BC spectrally and explicitly?

We appreciate any material to look at or any clue.


问:
Solving L1 minimization problems in Eigen
I have an $m\times n$ dense matrix $\mathbf{M}$ and wish to solve $\mathbf{M}\mathbf{x} = \mathbf{b}$ via any of the "L1" methods using Eigen. By this I mean I'm happy with using basis pursuit denoising, lasso, elasticnet, LARS, etc. But I would like it to be "canned", or at least somewhat easy using Eigen's primitives. The closest thing I can find are L2 methods, documented here.

Do there exist libraries building these solvers on top of Eigen? It appears from this slidedeck that a fundamental tool here is ADMM, yet I have not managed to find a canned implementation.


问:
Simulating Quantum Wave Function/Schrodinger Equation With A Time Varying Potential
I have solved the Time Independent Schrodinger Equation using the Numerov method and diagonalizing the Hamiltonian, in 1 - 3 dimensions. I suppose I could time-evolve it by multiplying every element in it by $e^{\frac{iEt}{\hbar}}$ to time evolve it (unless I cant). But I need an even more dynamical solution for what I need to do, and I need a time varying potential.

Is there some method where I can take an initial, known wave-function (like a free-particle solution) and evolve it in time while changing the potential? I have thought about using a Ritz solver, would that do the trick?

I would think the most practical solutions would be with the Schrodinger Equation, but if someone has a Quantum Field Theory solution that is also welcome.

Thank you!


问:
Partition of unity in FEM with bubble functions
In FEM with bubble functions, the field ($\boldsymbol{u}$) is approximated as a linear combination of the standard one ($\tilde{\boldsymbol{u}}$) plus the bubble field ($\boldsymbol{u}^b$). That is,
\begin{equation}
\boldsymbol{u} = \tilde{\boldsymbol{u}} + \boldsymbol{u}^b.
\end{equation}

Considering the addition of new basis functions and DOFs, I wonder about the partition of unity property.

Will it be violated? If yes, what are the potential consequences? If not, are the basis functions scaled accordingly?
Appreciate any insights and concerned references.


问:
function optimization using reinforcement learning
I'm regular with optimization process, but this time, I would like to optimize a multiobjective function with reinforcement learning, using python.

Of course I already have some testing functions (mathematical or more realistic), but I have no idea on how to proceed with RL.

I would like to compare results with more "traditional" approach.

If you have some tutorials or examples, it would be welcome.

Regards.


问:
Maximal "Convex Augmentation" of a Triangle in 2D Mesh
Consider a convex polygon in $\mathbb{R}^2$ with multiple convex holes in it and suppose that, for now, we have a 2D triangular mesh of the polygon, which is represented by $\mathcal{T} \equiv\{T_i\}_i$, a set of triangles in $\mathbb{R}^2$. My question is, given a triangle $T$ that lies at the exterior boundary of the polygon, is there a way to find a set $\{T_{i_j}\}_j \subset \mathcal{T}$ such that

$$
\left(\bigcup_j T_{i_j}\right)\cup T
$$

forms the largest convex polygon (without any holes) that contains $T$? Any help or hint will be greatly appreciated!

My Attempts:

I first wanted to solve the problem with a greedy algorithm. That is, check whether the union between $T$ and (not necessarily all) other triangles that share edges with $T$ is convex. If so, combine them and repeat this step with the resulting shape, and terminate when no more such triangles can be merged. But I find it hard to guarantee that the final resulting shape will be maximal.


问:
Are both centered and upwind numerical fluxes correct in DGTD?
In discontinuous Galerkin time domain (DGTD) method, a critical concept is the numerical flux that is used to link neighbouring elements. The numerical flux is however not unique. The popular choices include centered and upwind ones. I see in many articles and books describe their different numerical behaviours, but a problem has been hovering over my head for quite some time. Since distinct flux choices lead to DIFFERENT results, so rigorously speaking all these results CAN'T be correct. By "correct" I mean the solution faithfully satisfies the original partial differential equations. My question is: which flux choice is correct? (why do we need the wrong ones?) Or all the numerical fluxes are approximate? In the latter case, which is a better one?

Maybe my question is naive for people working with DGTD, the literature I read really did a good job to cause my confusion... Thanks for any comment!


问:
Matrix-free FEM references
I've seen that many people are using matrix-free fem codes in my community (mechanical engineering). I have to admit that I googled a bit and I didn't manage to find a good reference for the subject. Even a quick look to the books cited in this question Matrix free finite elements method for visualization in process tomography are not really introducing the matrix-free concept in the FEM context, rather they're giving it for granted. I'd like to know some references people have used as an introduction to the topic.</p>

Title: Matrix-free FEM references


问:
Looking for a specific version of the Quasi-Minimal Residual (QMR) method
I'm looking for an alternative formulation of Quasi-Minimal Residual (QMR) from Freund and Nachtigal (1994) based on a Lanczos process for complex valued matrices based on $A^H$ instead of $A^T$.

Specifically, I'm not quite sure how algorithm 7.1 would change if one where to use $A^H$ instead of $A^T$. Is there some article or book where I can find such a QMR formulation ?


问:
Losing memory on each call to scalapack solve function
I am solving a large system of equations using scalapack. Some systems that should run with no problem failed, apparently due to lack of memory in an mpi call. After investigating with execute_command_line("free -m"), I found that I'm losing large amounts of available RAM after every call to PZGETRS. By large, I mean O(GB), but it varies based on the matrix size and also the number of processes I include in the scalapack context. After some number of calls to PZGETRS, I do eventually run out of RAM, leading to the error messages.

I created a minimum working example:

```<code>program main
use mpi_f08
implicit none

    integer(kind=4) :: ierr, commW, commM, num_tasks, icontxt, nprow, npcol, n, &
        info, sub, trow, ii, my_row, my_col, tcol
    complex(kind=8), dimension(:,:), allocatable :: A, BX
    integer(kind=4), dimension(:), allocatable :: ipiv
    integer(kind=4), dimension(9) :: block_desc
    integer, external :: numroc   ! blacs routine
    double precision :: t1, t2, tf, ts

    ! Initialize MPI and establish self
    call mpi_init(ierr)
    commW=mpi_comm_world%mpi_val
    call mpi_comm_rank(mpi_comm_world, commM, ierr)
    call mpi_comm_size(mpi_comm_world, num_tasks, ierr)

    if (commM == 0) write(*,*) "Program start"
    if (commM == 0) call execute_command_line("free -m | head -2 | sed -e 1d")

    ! Initalize process array
    nprow = ceiling(sqrt(real(num_tasks,kind=4)))    ! Make grid squareish
    do while (mod(num_tasks, nprow) /= 0)
        nprow = nprow + 1
    end do
    npcol = num_tasks/nprow
    call blacs_get(0,0,icontxt)
    call blacs_gridinit(icontxt,'R',nprow,npcol)
    call blacs_gridinfo(icontxt, nprow, npcol, my_row, my_col)

    ! Set up arrays based on number of tasks and the small block size
    n = 6200
    sub = 40            ! size of block in block cyclic dist.
    trow = numroc(n,sub,my_row,0,nprow)
    tcol = numroc(n,sub,my_col,0,nprow)
    call descinit(block_desc, n, n, sub, sub, 0, 0, icontxt, trow, ierr)

    ! everyone allocate their data and fill with something
    allocate(A(trow, tcol), BX(trow, tcol), ipiv(n))
    A = cmplx(1.d0,1.d0,kind=8)
    BX = cmplx(2.d0,1.d0,kind=8)
    do ii = 1,min(trow,tcol)
        A(ii,ii) = cmplx(2.d0,2.d0,kind=8)
    end do

    ! call factor routine then solve
    if (commM == 0) write(*,*) "All of my variables allocated."
    if (commM == 0) call execute_command_line("free -m | head -2 | sed -e 1d")
    call mpi_barrier(mpi_comm_world)
    t1 = mpi_wtime()
    call pzgetrf(n,n,A,1,1,block_desc,ipiv,info)
    t2 = mpi_wtime()
    tf = t2 - t1
    if (commM == 0) write(*,*) "Post factor."
    if (commM == 0) call execute_command_line("free -m | head -2 | sed -e 1d")
    call mpi_barrier(mpi_comm_world)
    t1 = mpi_wtime()
    call pzgetrs('N',n,n,A,1,1,block_desc,ipiv,BX,1,1,block_desc,info)
    t2 = mpi_wtime()
    ts = t2 - t1
    if (commM == 0) write(*,*) "Post solve."
    if (commM == 0) call execute_command_line("free -m | head -2 | sed -e 1d")
    call mpi_barrier(mpi_comm_world)

    deallocate(A,BX,ipiv)
    call blacs_gridexit(icontxt)
    if (commM == 0) write(*,*) "Shutting down."
    if (commM == 0) call execute_command_line("free -m | head -2 | sed -e 1d")
    if (commM == 0) write(*,*) "Factor: ", tf, "Solve: ", ts
    call mpi_finalize(ierr)

end program main
</code>```

I'm linking with MKL according to the intel link line advisor and the cluster documentation that I'm using. Some lines from my Makefile, Makefile.inc...

```<code>LMPI                 := -L$(TACC_IMPI_LIB) -lmpi
IMPI                 := -I$(TACC_IMPI_INC)
LMKL                 := -L$(TACC_MKL_LIB) -lmkl_scalapack_lp64 -lmkl_blacs_intelmpi_lp64 -lpthread -lm -ldl
IMKL                 := -I$(TACC_MKL_INC)
LCOMBINE             := $(LMPI) $(LMKL)
ICOMBINE             := $(IMPI) $(IMKL)
...
FFLAGS := -O0 -g -debug -mkl=cluster
FC := $(TACC_IMPI_BIN)/mpiifort
FL := $(TACC_IMPI_BIN)/mpiifort
...
@echo Compiling $<
@$(FC) $(FFLAGS) $(ICOMBINE) -c $< -o $(BIN_DIR)/$@
@echo Linking to executable "'$(MAIN_EXE)'"
@$(FL) $(FFLAGS) $(MAIN_OBJ_FINAL) $(LCOMBINE) -o $(MAIN_EXE)
</code>```

(The full Makefiles are uploaded on my question for the intel forum, currently unsolved).

After lots of tests, I've found the memory loss depends on the software I compile my code with. I would prefer to use impi/intel-mkl given how quickly that combination solves my problems. However, that combo uses the most memory. Using OpenMPI + open source scalapack (v2.0.0), my memory loss is much smaller (for example, ~7300 MB with intel, ~200 MB with openmpi), but my code runs significantly slower.

Intel documentation mentions that some buffers are not deallocated until the program shuts down, for speed reasons, and says to either set the environmental variable MKL_DISABLE_FAST_MM=1 to turn that off, or have a line in the code: mm = mkl_disable_fast_mm() before calling mkl functions to disable this feature. Neither option has worked for me. I used the memory probe functions provided by intel to monitor maximum memory use, and got back a claim that 11 MB were used when I watch my available RAM drop by 7 GB via the free command.

My question is why am I'm seeing this behavior? Am I compiling wrong, or using scalapack wrong, or is this an actual bug? Does anyone have any suggestions to fix this?


问:
Over-specification of conjugate heat transfer coupling conditions
I am trying to implement steady state conjugate heat transfer using a monolithically coupled scheme. In this simulation, the computational domain is divided into fluid and solid subdomains. Over time, the fluid should heat up and transfer heat to the solid which is conducted within the solid. Heat is transferred by enforcing conservation of energy on the fluid-solid interface which requires that the heat flux is conserved via

$\mathbf{q}_f \cdot \mathbf{n} = \mathbf{q}_s \cdot \mathbf{n}$

where the subscripts $f$ and $s$ refer to quantities on the fluid and solid domain, respectively.

I am pursuing what the literature refers to perfect heat which assumes perfect contact (i.e. no gaps) between the solid and the fluid. One consequence of this is that we enforce the condition

$T_f = T_s$

on the fluid-solid interface. Both conditions are currently (and should ideally be) two-way coupled, meaning that information from both domains is employed to determine the correct state on the fluid-solid interface, as opposed to the state of one domain being prescribed from the state on the other domain.

My question is: are these coupling constraints too strict? My solver is able to satisfy one of these constraints in isolation, but never both. However, both constraints seem to be necessary.

A major difficulty seems to be that the time scale for heat conduction within the solid is much larger than heat convection in the fluid. The fluid heats up, supplying a large heat flux to the solid. However, the heat flux in the solid is very small due to the large time scale needed for conduction, ultimately reducing the solid heat flux and cooling the fluid. Requiring the temperatures to be equal would ideally increase the solid heat flux and prevent cooling of the fluid.

Any direction or papers are appreciated!


问:
Compact Finite Differences for the Heat Equation with Robin Boundary Conditions
I am trying to solve the Heat equation with Robin Boundary condition:

$$ u_t(x,t) = u_{xx}(x,t), \\ u(x,0) = g(x), \\ u(0,t) + u_x(0,t) = h_0(t), \\ u(1,t) + u_x(1,t) = h_1(t)$$

for $ 0\leq x\leq1$ and $t>0$, by using twice the composition of a fourth-order compact finite difference (CFD) approximation of the first derivative (taken from S. Lele (1991)) for the spatial discretization, i.e.

$$ M \; \partial_x \textbf{u} = A \;\textbf{u} \\ M \; \partial^2_x \textbf{u} = A \; \partial_x \textbf{u}, $$

where $M, A$ are two sparse matrices, corresponding to the coefficients obtained from the CFD scheme. This would simplify as follows

$$ \partial^2_x \textbf{u} = (M^{-1} A)^2 \;\textbf{u}. $$

As adviced in S. Lele (1991), I use an upwinding and downwinding scheme at the boundaries, so I avoid using ghost cells with the derivatives.

I am aware this could be done in an easier manner by using an approximation of the second derivative, as done in J. Malele et al. (2022), but my problem is in fact more complicated than the one shown here and requires me to use twice the composition of the first derivative.

And here, my question: How could I include the Robin boundary conditions?

I was advised to include the equations as additional lines on the last system, so that

$$ \hat{\textbf{u}} = \begin{bmatrix} u_x\Big|_1 \\ u_0 \\ u_1 \\ \vdots \\ u_N \\ u_x\Big|_{N-1} \end{bmatrix} $$

would be my new unknown, however, that makes a new matrix appear on the left hand side of the equation, i.e. left-multiplying $ \partial_x^2 \hat{\textbf{u}} $, which is a singular matrix. The fact that this is a singular matrix would then complicate the resolution of the ODE with another method time integrator.

Thank you for your time.

P.S.: The matrix $M$ is the matrix of coefficients of the derivatives. That means that if the schemes are of this form (as per S. Lele (1991)):

$$ \alpha f_{i-1}' + f_i + \alpha f_{i+1}' = \frac{a}{2h} \left( f_{i+1} - f_{i+1} \right)  \quad \forall i = 1,2,\dots, N-1$$

at the inner nodes and

$$ f_0' + \beta f_1' = \frac{1}{h}\left( bf_0 + cf_1 +df_2 + ef_3 \right) $$

at the left boundary, i.e. $i = 0$, and something similar for the right boundary, then the Matrix of coeffients will be:

$$ M := \begin{bmatrix} 1 & \beta & 0 & 0 & \dots \\ \alpha & 1 & \alpha & 0 & \ddots \\ 0 & \alpha & 1 & \alpha & \ddots \\ \vdots & \ddots & \ddots & \ddots & \ddots \end{bmatrix} $$


问:
Vehicle passenger assignment with capacity constraint

Problem Background

I'm trying to find a solution to the following passenger matching problem:

The network is represented by graph $G=(V,E)$. $V$ is the set of nodes/stations. $p_{ij}$ is the profit of traveling through an edge $(i,j)$. Let $N$ be the number of vehicles, all the vehicles have the same capacity $c$.

At each (discreet) time step, the passengers arrive at their origin station and need to be transported to their destination station. $d^t_{ij}$ is the demand and $f^t_{ij} \le d^t_{ij}$ is the passenger flow, i.e., the number of passengers that are traveling from $i$ to $j$ at time $t$ and are successfully matched to a vehicle. The unmatched passers will leave the system. $X^t_i$ is the total number of available vehicles at station $i$ at time $t$.

Objective

The objective is to maximize profit.

Problem formulation for single occupancy vehicle

I read a few resources and found the following formulation. However, this formulation assumes that the vehicles are single occupancy.

$$ \sum_{i,j \in V} f^t_{ij}p^t_{ij} $$
$$ 0 \le f^t_{ij} \le d^t_{ij} \quad i,j \in V \quad (1)$$
$$ \sum_{j\in V} f^t_{ij} \le X^t_i \quad i \in V \quad \quad (2)$$

The constraint $(1)$ ensures that the passenger flow doesn't exceed the demand and is non-negative.

The constraint $(2)$  ensures that the number of vehicles doesn't exceed the number of available vehicles $X^t_i$  at time $t$ at station $i$.

As long as the demand is integral, this formulation results in integral passenger flow.

Question

I'd like to extend the above single occupancy formulation to multi-passenger vehicles. I explored some resources online but the methods focused on individual vehicles instead of stations and were using ILP.

If anyone can point me to some related models and methods, or suggest how to solve this problem I would be grateful.

Update 1

I thought to add the capacity $c$ to the 2nd constraint as follow:

$$ \sum_{j\in V} f^t_{ij} \le X^t_{i}c \quad i \in V \quad \quad (2)$$

Because I have a limited knowledge of optimization problems, I'd appreciate any feedback regarding this.


问:
SUNDIALS using Python Interface tutorial
I am attempting to solve differential algebraic equations using SUNDIALS with a Python interface. I am having trouble learning to actually use sundials, and cannot seem to find any tutorials online that show some python examples. Would anyone have any links?</p>


问:
Why do BVP solvers' APIs only allow "unknown" parameters in the derivative and residual functions but not "known" parameters?
I recently needed to solve a second order boundary value problem and noticed that both scipy.integrate.solve_bvp and Matlab's bvp4c and bvp5c do not seem to support the passing of known parameters to the function odefun defining $\frac{dy}{dx}$ or the residual function which determines the boundary condtions bcfun.

They do however allow for unknown parameters to be simultaneously inferred while solving the problem. This seems like a strange API limitation, so I wonder if there is a reason for this.

I am used to working with solvers for initial value problems and optimizers for general minimization problems, and in both cases the APIs I'm familiar with allow defining the relevant functions such that they accept additional "known" parameters, e.g. scipy.integrate.solve_ivp allows defining
$\frac{dy}{dx} = f(x, y, a, b, c)$ and then passing the values for $a$, $b$ and $c$ as a tuple when solving the problem.


问:
Why can we remove the half-step velocity in velocity Verlet

Eliminating the half-step velocity, this algorithm may be shortend to

Why can we eliminate the half-step velocity and all the math behind the velocity Verlet to what Wikipedia shows?


问:
generating multivariate Gaussian samples with constraints
Using GRFs, one generates samples $x \sim N(m,K)$ with arbitrary mean $m$ and covariance matrix $K$ using a scalar Gaussian generator as

$x = m + Lu$,

where $L$ is a lower triangular matrix, such that $K =L^TL$ and $u \sim N(0,I)$.

I would like to generate samples, such that the first and last elements of x has a certain value.
My goal is to generate random right-hand sides for solving PDE, but such that the boundary conditions are fulfilled (1D).


问:
navier-stokes equations in semi-discrete form
Can someone point me to where the Navier-Stokes equations in 2D, both compressible and incompressible are written in semi-discrete form? I'm doing reduced order modelling and I need to write them down in the form

$$
\frac{d\mathbf{x}}{dt} = f(\mathbf{x},t,\mu)
$$

Where $\mu$ is some parameter and $\mathbf{x}$ is a high-dimensional vector containing velocities and pressures at all nodes in the mesh and $t$ is time.

I can write Burgers' equation in 2D in semi-discrete form; but with Navier-Stokes I'm confused about what to do with the pressure term. In incompressible Navier-Stokes we solve a linear system for the pressure and I'm not sure about how to include that in the form I wrote above.

I'm confused because in incompressible Navier-Stokes the time-derivative appears only in the velocity and, in the only algorithm I'm familiar with, the pressure projection algorithm there is a linear equation for the pressure. These steps mean that I can't see how to write the Navier-Stokes equations in the form I need.


问:
solve a coupled PDE system with some discontinuity by a mixed FEM
$$
\begin{cases}
u_t=f(u,v,\nabla{v}),\\
\Delta{[g(u)v]}=0.
\end{cases}$$
I want to solve the above PDE system by a mixed FEM, that is, $u_t=f(u,v,\nabla{v})$ by discontinuous Galerkin (DG), where $f$ doesn't depend on $\nabla{u}$, and $\Delta{[g(u)v]}=0$ by continuous Galerkin (CG). The key is that in the true solution of such a PDE system, $u$ and $\nabla{v}$ exhibit some spatial discontinuity/shock but $v$ is spatially continuous. For this reason, I'd like to solve $u_t=f(u,v,\nabla{v})$ by DG to retain the spatial discontinuity of both $u$ and $\nabla{v}$.

As is known to all, DG was first used to solve $u_t+u_x=0$. Since $u_t=f(u,v,\nabla{v})$ isn’t the standard form of $u_t+u_x=0$, how can I implement the DG scheme here given the complexity of the true solution?


问:
What is the boundary condition of $\varepsilon $ for $k-\varepsilon $ turbulence model?
I'm working on the numerical solution of systems of equations, you can use this link to access it, the system of equations is
$\begin{aligned} \partial_t \theta+u \nabla \theta-\nabla \cdot\left(\kappa_T^m \nabla \theta\right) & =0 \\ \partial_t u+u \nabla u-\nabla \cdot\left(\mu_T \nabla u\right)+\nabla p+e\left(\theta-\theta_0\right) \vec{e}_2 & =0 \\ \nabla \cdot u & =0 \\ \mu_T & =c_\mu \frac{k^2}{\epsilon} \\ \kappa_T & =\kappa \mu_T \\ \partial_t k+u \nabla k+\epsilon-\nabla \cdot\left(\mu_T \nabla k\right) & =\frac{\mu_T}{2}\left|\nabla u+\nabla u^T\right|^2 \\ \partial_t \epsilon+u \nabla \epsilon+c_2 \frac{\epsilon^2}{k}-\frac{c_\epsilon}{c_\mu} \nabla \cdot\left(\mu_T \nabla \epsilon\right) & =\frac{c_1}{2} k\left|\nabla u+\nabla u^T\right|^2\end{aligned}$
when dealing with the last equation, I need to deal with the normal derivative of $\varepsilon$ on a boundary of wall, could you please tell me what is the normal derivative of $\varepsilon$ on the impermeable wall? thank you.


问:
Exponential Integrator to solve PDE with Stiff term
I wish to solve an equation like the following,
$$\frac{\partial f}{\partial t}+\frac{\partial}{\partial x}\left(A(x)f\right)=S(x,t)f$$
where $A(x,f)f$ and $S(x,t)f$ are the advection and the source terms.
When I naively try to incorporate the source term while discretizing the advection part via an upwind scheme, for example, the solution starts to generate oscillations.
The discretization is described in the following,
$$\frac{f^{n+1}_{i}-f^{n}_{i}}{dt}=-\frac{\mathcal{F}_{i+\frac{1}{2}}^{n}-\mathcal{F}_{i-\frac{1}{2}}^{n}}{x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}}+S\left(x_{i},t^{n}\right)f^{n}_{i}$$
with $$\mathcal{F}_{i+\frac{1}{2}}^{n}=\begin{cases}
      A\left(x_{i+\frac{1}{2}}\right)f^{n}_{i} & \text{if $A\left(x_{i+\frac{1}{2}}\right)>0$}\\
      A\left(x_{i+\frac{1}{2}}\right)f^{n}_{i+1} & \text{if $A\left(x_{i+\frac{1}{2}}\right)<0$}
    \end{cases}. $$
The same discretization is giving the right output when $S=0$.
I found out the problem is related to the time dependency of $S$, the source is sometimes getting very large which is causing the problem.
I found out that exponential integrator can be used to mitigate such problems but I am not quite sure how to implement it for the current problem.
The literature that I find starts with the following equation $$\frac{df}{dt}=\mathcal{K}f+\mathcal{G}(t,f)\\
\implies f^{n+1}=e^{\mathcal{K}dt}+\int_{0}^{dt}e^{\mathcal{K}\left(dt-\tau\right)}G\left(t^{n}+\tau,f^{n+\tau}\right)\,d\tau$$
where $\mathcal{K}$ is the term responcible for stiffness and $\mathcal{G}$ is some function.
I understand for the current problem $\mathcal{K}=S(x,t)$ but I do not know how to connect $\mathcal{G}$ with the advection part and how to compute the above integral.
Any help in this regard will be highly appreciated.


问:
books/paper recommendation on computational thermal-turbulence by using FEM
I have just learned basic FEM for 2D N-S euqation, now my teacher let me to do the following problem, the document of this problem is in large fluid problem, the system of equations is listed in that page, I don't even know what is the boundary condition, could you please recommend some textbooks or papers for solving the following system of equations by using FEM?thank you very much.$$
\begin{aligned}
\partial_t \theta+u \nabla \theta-\nabla \cdot\left(\kappa_T^m \nabla \theta\right) & =0 \\
\partial_t u+u \nabla u-\nabla \cdot\left(\mu_T \nabla u\right)+\nabla p+e\left(\theta-\theta_0\right) \vec{e}_2 & =0 \\
\nabla \cdot u & =0 \\
\mu_T & =c_\mu \frac{k^2}{\epsilon} \\
\kappa_T & =\kappa \mu_T \\
\partial_t k+u \nabla k+\epsilon-\nabla \cdot\left(\mu_T \nabla k\right) & =\frac{\mu_T}{2}\left|\nabla u+\nabla u^T\right|^2 \\
\partial_t \epsilon+u \nabla \epsilon+c_2 \frac{\epsilon^2}{k}-\frac{c_\epsilon}{c_\mu} \nabla \cdot\left(\mu_T \nabla \epsilon\right) & =\frac{c_1}{2} k\left|\nabla u+\nabla u^T\right|^2
\end{aligned}
$$
I want to know the boundary conditions, choose of element, method of transform the nonlinear term to linear term,could you please recommend some books or papers, thank you


问:
Fortran - Lid-Driven Cavity Boundary Conditions Error when using SIMPLE method
I am studying Numerical Methods for incompressible flows. part of the tasks is to model the lid driven cavity problem in 2D using the SIMPLE method.

I have been provided with Fortran code that is working solved, however is for a channel flow problem. I was informed that by modifying the boundary conditions I can convert it from channel flow to cavity flow. (by manipulating which walls as velocity etc etc).

I have created a SIMPLE Solver for the the Lid Driven Cavity Problem in MATLAB. It is known working and Validated against the Ghia et al.(1982) However when I modify the fortran code to what I believe to be the correct boundary conditions my results.dat file contains NaN for all values of U and V for all X and Y locations.

here is the original Fortran code (in its channel form).

```program SIMPLE_TDMA_2D
! variables
! dimensions:
      integer :: iNxCells, iNyCells, iNxNodes, iNyNodes, iNxUNodes, iNyUnodes, &
                 iNxVNodes,iNyVNodes,iNxPNodes,iNyPNodes, iNxMax, iNyMax
!iterations
      integer :: iMaxIter, i,j,n
!convergence limit
      real :: rConvergence
!residuals
      real :: rURes, rVRes, rPRes, rMaxRes, rCRes
!underrelaxation
      real :: rAlphaP, rAlphaU, rAlphaV, rAlphaC
!Dimensionless numbers
      real :: rReynolds, rPeclet, rInjStart,rInjEnd
!massflow
      real :: rTotFlow, rMRatio, rCInFlow, rCOutFlow
!grid
      real :: rHeight, rWidth, rDX, rDY
! current auxiliary variables
      real :: rPCur, rUCur,rVCur,rCCur
! storage for variables
      real, allocatable, dimension (:,:) :: rU,rV,rP,rC,rDU,rDV,rDP

! read parameters from file

      open (unit=1, file="parameters.par",status="old")
      read (1,*) rHeight
      read (1,*) rWidth
      read (1,*) rReynolds
      read (1,*) rPeclet
      read (1,*) rInjStart
      read (1,*) rInjEnd
      read (1,*) iNxCells
      read (1,*) iNyCells
      read (1,*) iMaxIter
      read (1,*) rConvergence
      read (1,*) rAlphaP
      read (1,*) rAlphaU
      read (1,*) rAlphaV
      read (1,*) rAlphaC
      close(1)
! grid cells
      rDx=rWidth/iNxCells
      rDy=rHeight/iNyCells
! grid nodes
      iNxNodes=iNxCells+1
      iNyNodes=iNyCells+1
! nodes for P - control volumes + 2 for BCs
      iNxPNodes=iNxCells+2
      iNyPNodes=iNyCells+2
! nodes for U - control volumes + 1 in X and 2 in Y
      iNxUNodes=iNxCells+1
      iNyUNodes=iNyCells+2
! nodes for V - control volumes + 2 in X and 1 in Y
      iNxVNodes=iNxCells+2
      iNyVNodes=iNyCells+1
!
      iNxMax=iNxCells+2
      iNyMax=iNyCells+2
! allocate memory
! scalar - same control volume as P
      allocate(rU(iNxMax,iNyMax),rV(iNxMax,iNyMax),rP(iNxMax,iNyMax),rC(iNxMax,iNyMax))
      allocate(rDU(iNxMax,iNyMax),rDV(iNxMax,iNyMax),rDP(iNxMax,iNyMax))
! Initialise
      rP=0.
      rU=0.
      rV=0.
      rC=0.
      rMaxRes=10.e+3
      n=1
! start time loop
      do while(n.le.iMaxIter.and.rMaxRes.gt.rConvergence)
! compute u-momentum equation
       call UMomentum(iNxUNodes, iNyUNodes,iNxMax,iNyMax,rDx,rDy,rU,rV,rP,rDU,rAlphaU,rURes,rReynolds)
! compute v-momentum equation
       call VMomentum(iNxVNodes, iNyVNodes,iNxMax,iNyMax,rDx,rDy,rU,rV,rP,rDV,rAlphaV,rVRes,rReynolds)
! correct mass flux
       rTotFlow=0.
       do j = 2 , iNyUNodes-1
          rTotFlow=rTotFlow+rDY*rU(iNxUNodes,j)
       end do
       if(rTotFlow.gt.0.0) then
            rMRatio=rHeight/rTotFlow
         else
            rMRatio=1.0
         end if
       do j = 2 , iNyUNodes-1
          rU(iNxUNodes,j)=rU(iNxUNodes,j)*rMRatio
       end do
! check scalar flux
       rCOutFlow=0.
       do j = 2 , iNyUNodes-1
          rCOutFlow=rCOutFlow+rDY*rC(iNxUNodes,j)
       end do
! check scalar flux
       rCInFlow=0.
       do j = 2 , iNyUNodes-1
          rCInFlow=rCInFlow+rDY*rC(2,j)
       end do
!       print *, "Scalar: ", rCInFlow,rCOutFlow
! compute pressure correction
       call PCorrection(iNxPNodes, iNyPNodes,iNxMax,iNyMax,rDx,rDy,rU,rV,rDU,rDV,rDP,rPRes)
! correct velocities
       do i=2,iNxUNodes-1
        do j=2,iNyUNodes-1
          rU(i,j)=rU(i,j)-rDU(i,j)*(rDP(i+1,j)-rDP(i,j))
        end do
       end do
       do i=2,iNxVNodes-1
        do j=2,iNyVNodes-1
          rV(i,j)=rV(i,j)-rDV(i,j)*(rDP(i,j+1)-rDP(i,j))
        end do
       end do
! correct pressure
       rP=rP+rAlphaP*rDP
! scale pressure
       rP=rP-rP(2,2)
! Solve the scalar equation for corrected flow field:
       call Scalar(iNxPNodes, iNyPNodes,iNxMax,iNyMax,rDx,rDy,rU,rV,rC,rAlphaC,rCRes,rPeclet,rInjStart,rInjEnd)
! maximum residual
       rMaxRes=max(rURes,rVRes,rPRes,rCRes)
       write(*,"(2X,I4,2X,5(2X,E10.3))") n,rURes,rVRes,rPRes,rCRes,rP(2,2)
! advance iterations
       n=n+1
      end do
! check convergence3
      if(rMaxRes.gt.rConvergence) then
       write(*,*) "Stat: Convergence not reached. Exiting..."
      else
       write(*,*) "Stat: Converged."
      end if
!output tecplot file:
      open (unit=1,file="resultn.dat",status="replace",action="write")
      write(1,*) 'Variables="X","Y","P","U","V", "Scalar"'
      write(1,'("ZONE DATAPACKING=POINT, I=",I4,", J=",I4)') iNxNodes,iNyNodes
      do j=1,iNyNodes
       do i=1,iNxNodes
         rPCur=0.25*(rP(i,j)+rP(i+1,j)+rP(i,j+1)+rP(i+1,j+1))
         rCCur=0.25*(rC(i,j)+rC(i+1,j)+rC(i,j+1)+rC(i+1,j+1))
         rUCur=0.5 *(rU(i,j)+rU(i,j+1))
         rVCur=0.5 *(rV(i,j)+rV(i+1,j))
         write(1,"(8E15.7)") (i-1)*rDX,(j-1)*rDY,rPCur,rUCur,rVCur,rCCur
       end do
      end do
      close(1)
end program SIMPLE_TDMA_2D

subroutine UMomentum(iNxUNodes, iNyUNodes,iNxMax,iNyMax,rDx,rDy,rU,rV,rP,rDU,rAlphaU,rURes,rReynolds)
!arguments
    integer, intent(in) :: iNxUNodes, iNyUNodes,iNxMax,iNyMax
    real, intent(in)    :: rDx,rDy, rAlphaU,rReynolds
    real, intent(inout) :: rURes
    real, dimension (iNxMax, iNyMax), intent (in) :: rV,rP
    real, dimension (iNxMax, iNyMax), intent (inout) :: rU,rDU
! local variables
    real ::  de, dw, dn, ds, ge, gw, gn, gs,rTemp
    real, dimension(iNxMax,iNyMax) ::aw,ae,as,an,ap,su
    integer :: i,j

! set boundary conditions for coefficients in the equations
! south
    aw(:,1)=0.
    ae(:,1)=0.
    as(:,1)=0.
    an(:,1)=0.
    ap(:,1)=1.
    su(:,1)=0.
! west
    aw(1,:)=0.
    ae(1,:)=0.
    as(1,:)=0.
    an(1,:)=0.
    ap(1,:)=1.
    su(1,:)=1.
! north
    aw(:,iNyUNodes)=0.
    ae(:,iNyUNodes)=0.
    as(:,iNyUNodes)=0.
    an(:,iNyUNodes)=0.
    ap(:,iNyUNodes)=1.
    su(:,iNyUNodes)=0.
! east
    aw(iNxUNodes,:)=1.
    ae(iNxUNodes,:)=0.
    as(iNxUNodes,:)=0.
    an(iNxUNodes,:)=0.
    ap(iNxUNodes,:)=1.
    su(iNxUNodes,:)=0.
! set the coefficeints in the interior staggered cells

      do i = 2 , iNxUNodes-1
         do J = 2 , iNyUNodes-1
            dw = rDY/(rDX*rReynolds)
            gw = rDY*(rU(i-1, j) + rU(i, j)) / 2.
            de = rDY/(rDX*rReynolds)
            ge = rDY*(rU(i, j) + rU(i+1, j)) / 2.
            if ( j .eq. 2 ) then
               ds = 2.*rDX/(rDY*rReynolds)
               gs = 0.0
            else
               ds = rDX/(rDY*rReynolds)
               gs = rDX*(rV(i, j-1) + rV(i+1, j-1)) / 2.
            end if
            if ( j .eq. (iNyUNodes - 1)) then
               dn = 2.*rDX/(rDY*rReynolds)
               gn = 0.0
            else
               dn = rDX/(rDY*rReynolds)
               gn = rDX*(rV(i, j) + rV(i+1, j)) / 2.
            end if
            aw(i, j) = dw + max( 0.0 ,  gw)
            ae(i, j) = de + max( 0.0 , -ge)
            as(i, j) = ds + max( 0.0 ,  gs)
            an(i, j) = dn + max( 0.0 , -gn)
            ap(i, j) = dw+de+ds+dn+max( 0.0 ,  -gw) + max( 0.0 , ge) +max( 0.0 ,  -gs)+max( 0.0 ,  gn)
            su(i, j) = -(rP(i+1, j) - rP(i, j))*rDY
         end do
      end do

! evaluate the residual in the interior
      rURes=0.
      do i = 2 , iNxUNodes - 1
         do j = 2 , iNyUNodes - 1
           rTemp = abs( ap(i, j)*rU(i, j) &
                      - ae(i, j)*rU(i+1, j) - aw(i, j)*rU(i-1, j) &
                      - an(i, j)*rU(i, j+1) - as(i, j)*rU(i, j-1) &
                      - su(i, j))
           if(rTemp.ge.rURes) rURes=rTemp
         end do
      end do
! apply underrelaxation to interior cells
      ap(2:iNxUNodes-1,2:iNyUNodes-1)=ap(2:iNxUNodes-1,2:iNyUNodes-1)/rAlphaU
      su(2:iNxUNodes-1,2:iNyUNodes-1)=su(2:iNxUNodes-1,2:iNyUNodes-1)+ &
         (1.-rAlphaU)*ap(2:iNxUNodes-1,2:iNyUNodes-1)*rU(2:iNxUNodes-1,2:iNyUNodes-1)
! store dU for pressure correction:
      rDU(2:iNxUNodes-1,2:iNyUNodes-1)=rDY/ap(2:iNxUNodes-1,2:iNyUNodes-1)
! solve the tridiagonal systems
!  for each i-direction line
      call tridag_i(as,aw,ap,ae,an,su,rU,iNxUNodes,iNyUNodes,iNxMax,iNyMax)
!  for each j-direction line
      call tridag_j(aw,as,ap,an,ae,su,rU,iNyUNodes,iNxUNodes,iNxMax,iNyMax)

      return
end subroutine UMomentum

subroutine VMomentum(iNxVNodes, iNyVNodes,iNxMax,iNyMax,rDx,rDy,rU,rV,rP,rDV,rAlphaV,rVRes,rReynolds)
!arguments
    integer, intent(in) :: iNxVNodes, iNyVNodes,iNxMax,iNyMax
    real, intent(in)    :: rDx,rDy, rAlphaV,rReynolds
    real, intent(inout) :: rVRes
    real, dimension (iNxMax, iNyMax), intent (in) :: rU,rP
    real, dimension (iNxMax, iNyMax), intent (inout) :: rV,rDV
! local variables
    real ::  de, dw, dn, ds, ge, gw, gn, gs, rTemp
    real, dimension(iNxMax,iNyMax) ::aw,ae,as,an,ap,su
    integer :: i,j

! set boundary conditions for coefficients in the equations
! south
    aw(:,1)=0.
    ae(:,1)=0.
    as(:,1)=0.
    an(:,1)=0.
    ap(:,1)=1.
    su(:,1)=0.
! west
    aw(1,:)=0.
    ae(1,:)=0.
    as(1,:)=0.
    an(1,:)=0.
    ap(1,:)=1.
    su(1,:)=0.
! north
    aw(:,iNyVNodes)=0.
    ae(:,iNyVNodes)=0.
    as(:,iNyVNodes)=0.
    an(:,iNyVNodes)=0.
    ap(:,iNyVNodes)=1.
    su(:,iNyVNodes)=0.
! east
    aw(iNxVNodes,:)=1.
    ae(iNxVNodes,:)=0.
    as(iNxVNodes,:)=0.
    an(iNxVNodes,:)=0.
    ap(iNxVNodes,:)=1.
    su(iNxVNodes,:)=0.
! set the coefficeints in the interior staggered cells

      do i = 2 , iNxVNodes-1
         do j = 2 , iNyVNodes-1
            if ( i .eq. 2 ) then
             dw = 2.*rDY/(rDX*rReynolds)
             gw = rDY*(rU(i-1, j) + rU(i-1, j+1)) / 2.
            else
             dw = rDY/(rDX*rReynolds)
             gw = rDY*(rU(i-1, j) + rU(i-1, j+1)) / 2.
            end if
            if ( i .eq. (iNxVNodes - 1)) then
             de = 2.*rDY/(rDX*rReynolds)
             ge = rDY*(rU(i, j) + rU(i, j+1)) / 2.
        else
             de = rDY/(rDX*rReynolds)
             ge = rDY*(rU(i, j) + rU(i, j+1)) / 2.
        end if
            ds = rDX/(rDY*rReynolds)
            gs = rDX*(rV(i, j-1) + rV(i, j)) / 2.
            dn = rDX/(rDY*rReynolds)
            gn = rDX*(rV(i, j) + rV(i, j+1)) / 2.

            aw(i, j) = dw + max( 0.0 ,  gw)
            ae(i, j) = de + max( 0.0 , -ge)
            as(i, j) = ds + max( 0.0 ,  gs)
            an(i, j) = dn + max( 0.0 , -gn)
            ap(i, j) = dw+de+ds+dn+max( 0.0 ,  -gw) + max( 0.0 , ge) +max( 0.0 ,  -gs)+max( 0.0 ,  gn)
            su(i, j) = -(rP(i, j+1) - rP(i, J))*rDX
         end do
      end do
! evaluate the residual in the interior
      rVRes=0.
      do i = 2 , iNxVNodes - 1
         do j = 2 , iNyVNodes - 1
           rTemp=abs( ap(i, j)*rV(i, j) &
                    - ae(i, j)*rV(i+1, j) - aw(i, j)*rV(i-1, j) &
                    - an(i, j)*rV(i, j+1) - as(i, j)*rV(i, j-1) &
                    - su(i, j))
           if(rTemp.ge.rVRes) rVRes=rTemp
         end do
      end do

!output tecplot file:
! apply underrelaxation to interior cells
      ap(2:iNxVNodes-1,2:iNyVNodes-1)=ap(2:iNxVNodes-1,2:iNyVNodes-1)/rAlphaV
      su(2:iNxVNodes-1,2:iNyVNodes-1)=su(2:iNxVNodes-1,2:iNyVNodes-1)+ &
        (1.-rAlphaV)*ap(2:iNxVNodes-1,2:iNyVNodes-1)*rV(2:iNxVNodes-1,2:iNyVNodes-1)
! store dV for pressure correction:
      rDV(2:iNxVNodes-1,2:iNyVNodes-1)=rDX/ap(2:iNxVNodes-1,2:iNyVNodes-1)
! solve the tridiagonal systems
!  for each i-direction line
      call tridag_i(as,aw,ap,ae,an,su,rV,iNxVNodes,iNyVNodes,iNxMax,iNyMax)
!  for each j-direction line
      call tridag_j(aw,as,ap,an,ae,su,rV,iNyVNodes,iNxVNodes,iNxMax,iNyMax)
      return
end subroutine VMomentum

subroutine PCorrection(iNxPNodes,iNyPNodes,iNxMax,iNyMax,rDx,rDy,rU,rV,rDU,rDV,rDP,rPRes)
!arguments
    integer, intent(in) :: iNxPNodes, iNyPNodes,iNxMax,iNyMax
    real, intent(in)    :: rDx,rDy
    real, intent(inout) :: rPRes
    real, dimension (iNxMax, iNyMax), intent (in) :: rU,rV,rDU,rDV
    real, dimension (iNxMax, iNyMax), intent (inout) :: rDP
! local variables
    real, dimension(iNxMax,iNyMax) ::aw,ae,as,an,ap,su
    integer :: i,j
    real :: rTemp

! initalise pressure correction
    rDP=0.
! set boundary conditions for coefficients in the equations
    aw=0.
    ae=0.
    as=0.
    an=0.
    ap=0.
    su=0.
! south
    aw(2:iNxPNodes-1,1)=0.
    ae(2:iNxPNodes-1,1)=0.
    as(2:iNxPNodes-1,1)=0.
    an(2:iNxPNodes-1,1)=1.
    ap(2:iNxPNodes-1,1)=1.
    su(2:iNxPNodes-1,1)=0.
! west
    aw(1,2:iNyPNodes-1)=0.
    ae(1,2:iNyPNodes-1)=1.
    as(1,2:iNyPNodes-1)=0.
    an(1,2:iNyPNodes-1)=0.
    ap(1,2:iNyPNodes-1)=1.
    su(1,2:iNyPNodes-1)=0.
! north
    aw(2:iNxPNodes-1,iNyPNodes)=0.
    ae(2:iNxPNodes-1,iNyPNodes)=0.
    as(2:iNxPNodes-1,iNyPNodes)=1.
    an(2:iNxPNodes-1,iNyPNodes)=0.
    ap(2:iNxPNodes-1,iNyPNodes)=1.
    su(2:iNxPNodes-1,iNyPNodes)=0.
! east
    aw(iNxPNodes,2:iNyPNodes-1)=1.
    ae(iNxPNodes,2:iNyPNodes-1)=0.
    as(iNxPNodes,2:iNyPNodes-1)=0.
    an(iNxPNodes,2:iNyPNodes-1)=0.
    ap(iNxPNodes,2:iNyPNodes-1)=1.
    su(iNxPNodes,2:iNyPNodes-1)=0.
! set the coefficeints in the interior staggered cells

      do j = 2 , iNyPNodes-1
       do i = 2 , iNxPNodes-1
           aw(i, j) = rDY*rDU(i-1, j)
           ae(i, j) = rDY*rDU(i, j)
           as(i, j) = rDX*rDV(i, j-1)
           an(i, j) = rDX*rDV(i, j)
           ap(i,j)  = aw(i, j)+ae(i, j)+as(i, j)+an(i, j)
           su(i, j) = rDY*(rU(i-1, j)-rU(i, j))+ &
                      rDX*(rV(i, j-1)-rV(i, j))
         end do
      end do
! sweep through the domain repeatedly and solve the algebraic equations in each line
      do i=1,10
        call tridag_i(as,aw,ap,ae,an,su,rDP,iNxPNodes,iNyPNodes,iNxMax,iNyMax)
        call tridag_j(aw,as,ap,an,ae,su,rDP,iNyPNodes,iNxPNodes,iNxMax,iNyMax)
      end do
! compute the pressure residual
      rPRes=0.
      do j=2,iNyPNodes-1
       do i=2,iNxPNodes-1
        rTemp=abs(su(i,j))
        if(rTemp.ge.rPRes) rPRes=rTemp

       end do
      end do
      return
end subroutine PCorrection

subroutine Scalar(iNxPNodes, iNyPNodes,iNxMax,iNyMax,rDx,rDy,rU,rV,rC,rAlphaC,rCRes,rPeclet,rInjStart,rInjEnd)
!arguments
    integer, intent(in) :: iNxPNodes, iNyPNodes,iNxMax,iNyMax
    real, intent(in)    :: rDx,rDy, rAlphaC,rPeclet,rInjStart,rInjEnd
    real, intent(inout) :: rCRes
    real, dimension (iNxMax, iNyMax), intent (in) :: rU,rV
    real, dimension (iNxMax, iNyMax), intent (inout) :: rC
! local variables
    real ::  de, dw, dn, ds, ge, gw, gn, gs, rTemp
    real, dimension(iNxMax,iNyMax) ::aw,ae,as,an,ap,su
    integer :: i,j

! set boundary conditions for coefficients in the equations
! south
    aw(:,1)=0.
    ae(:,1)=0.
    as(:,1)=0.
    an(:,1)=1.
    ap(:,1)=1.
    su(:,1)=0.
! west
    aw(1,:)=0.
    ae(1,:)=0.
    as(1,:)=0.
    an(1,:)=0.
    do j=2,iNyPNodes-1
     rTemp=(j-2)*rDY+0.5*rDY
     if (rTemp.ge.rInjStart.and.rTemp.le.rInjEnd) then
      ap(1,j)=1.
      su(1,j)=1.
     else
      ap(1,j)=1.
      su(1,j)=0.
     end if
    end do
! north
    aw(:,iNyPNodes)=0.
    ae(:,iNyPNodes)=0.
    as(:,iNyPNodes)=1.
    an(:,iNyPNodes)=0.
    ap(:,iNyPNodes)=1.
    su(:,iNyPNodes)=0.
! east
    aw(iNxPNodes,:)=1.
    ae(iNxPNodes,:)=0.
    as(iNxPNodes,:)=0.
    an(iNxPNodes,:)=0.
    ap(iNxPNodes,:)=1.
    su(iNxPNodes,:)=0.
! set the coefficeints in the interior staggered cells

      do i = 2 , iNxPNodes-1
         do j = 2 , iNyPNodes-1
          dw = rDY/(rDX*rPeclet)
          gw = rDY*rU(i-1, j)
          de = rDY/(rDX*rPeclet)
          ge = rDY*rU(i, j)
          if (j.eq.2) then
            ds = 0.0
            gs = 0.0
          else
            ds = rDX/(rDY*rPeclet)
            gs = rDX*rV(i, j-1)
          end if
          if (j.eq.(iNyPNodes-1)) then
           dn = 0.0
           gn = 0.0
          else
           dn = rDX/(rDY*rPeclet)
           gn = rDX*rV(i, j)
          end if
           aw(i, j) = dw + max( 0.0 ,  gw)
       ae(i, j) = de + max( 0.0 , -ge)
       as(i, j) = ds + max( 0.0 ,  gs)
       an(i, j) = dn + max( 0.0 , -gn)
       ap(i, j) =  dw+de+ds+dn+max( 0.0 ,  -gw) + max( 0.0 , ge) +max( 0.0 ,  -gs)+max( 0.0 ,  gn)
       su(i, j) = 0.0
         end do
      end do
! evaluate the residual in the interior
      rCRes=0.
      do i = 2 , iNxPNodes - 1
         do j = 2 , iNyPNodes - 1
           rTemp=abs( ap(i, j)*rC(i, j) &
                    - ae(i, j)*rC(i+1, j) - aw(i, j)*rC(i-1, j) &
                    - an(i, j)*rC(i, j+1) - as(i, j)*rC(i, j-1) &
                    - su(i, j))
           if(rTemp.ge.rCRes) rCRes=rTemp
         end do
      end do

!output tecplot file:
! apply underrelaxation to interior cells
      ap(2:iNxPNodes-1,2:iNyPNodes-1)=ap(2:iNxPNodes-1,2:iNyPNodes-1)/rAlphaC
      su(2:iNxPNodes-1,2:iNyPNodes-1)=su(2:iNxPNodes-1,2:iNyPNodes-1)+ &
        (1.-rAlphaC)*ap(2:iNxPNodes-1,2:iNyPNodes-1)*rC(2:iNxPNodes-1,2:iNyPNodes-1)
! solve the tridiagonal systems
!  for each i-direction line
      call tridag_i(as,aw,ap,ae,an,su,rC,iNxPNodes,iNyPNodes,iNxMax,iNyMax)
!  for each j-direction linE
      call tridag_j(aw,as,ap,an,ae,su,rC,iNyPNodes,iNxPNodes,iNxMax,iNyMax)
      return
end subroutine Scalar

!-----------------------------------------------------------------------
!     Tri-diagonal matrix solver (based on Numerical Recipes)
!
!     Solving i-lines, looping over domain in j-direction
!
!     Note signs of Coefficients
! note - r - righthand side
!-----------------------------------------------------------------------
      subroutine tridag_i(d, a, b, c, e, r, u, n, l,iNxMax,iNyMax)
!-----------------------------------------------------------------------
!     Argument type declarations
!-----------------------------------------------------------------------
      integer, intent(in) ::  n, l, iNxMax, iNyMax
      real, dimension(iNxMax,iNyMax) ::     a, b, c,d,e,r,u
!-----------------------------------------------------------------------
!     Local variable type declarations
!-----------------------------------------------------------------------
      integer j, m
      real    bet,gam(n),rhs
!-----------------------------------------------------------------------
!     Loop in j-direction
!-----------------------------------------------------------------------
      do j = 2 , l-1
!-----------------------------------------------------------------------
!     Forward recurrence
!-----------------------------------------------------------------------
         bet  = b(1, j)
         u(1, j) = r(1, j)/bet
         do m = 2 , n
            gam(m) = -c(m-1, j) / bet
            bet    = b(m, j) + a(m, j)*gam(m)
            rhs    = r(m,j) + d(m,j)*u(m,j-1) + e(m,j)*u(m,j+1)
            u(m, j)= (rhs + a(m, j)*u(m-1, j))/bet
         end do
!-----------------------------------------------------------------------
!     Backward recurrence
!-----------------------------------------------------------------------
         do m = n-1, 1, -1
            u(m, j) = u(m, j) - gam(m+1)*u(m+1, j)
         end do
!-----------------------------------------------------------------------
!     Next j line
!-----------------------------------------------------------------------
      end do
!-----------------------------------------------------------------------
!     finished so return
!-----------------------------------------------------------------------
      return
      end

!-----------------------------------------------------------------------
!     Tri-diagonal matrix solver (based on Numerical Recipes)
!
!     Solving j-lines, looping over domain in i-direction
!
!     Note signs of Coefficients
! note - r - righthand side
!-----------------------------------------------------------------------
      subroutine tridag_j(d, a, b, c, e, r, u, n, l,iNxMax,iNyMax)
!-----------------------------------------------------------------------
!     Argument type declarations
!-----------------------------------------------------------------------
      integer, intent(in) ::  l, n,iNxMax,iNyMax
      real, dimension (iNxMax,iNyMax) ::  a,b,c,d,e,r,u

!-----------------------------------------------------------------------
!     Local variable type declarations
!-----------------------------------------------------------------------
      integer i, m
      real    bet,gam(n),rhs
!-----------------------------------------------------------------------
!     Loop in i-diection
!-----------------------------------------------------------------------
      do i = 2, l-1
!
!-----------------------------------------------------------------------
!     Forward recurrence
!-----------------------------------------------------------------------
         bet  = b(i, 1)
         u(i, 1) = r(i, 1)/bet
         do m = 2, n
            gam(m) = -c(i, m-1) / bet
            bet    = b(i, m) + a(i, m)*gam(m)
            rhs    = r(i,m) + d(i,m)*u(i-1,m) + e(i,m)*u(i+1,m)
            u(i, m)= (rhs + a(i, m)*u(i, m-1))/bet
         end do
!-----------------------------------------------------------------------
!     Backward recurrence
!-----------------------------------------------------------------------
         do m = n-1, 1, -1
            u(i, m) = u(i, m) - gam(m+1)*u(i, m+1)
         end do
!-----------------------------------------------------------------------
!     Next i line
!-----------------------------------------------------------------------
      end do
!-----------------------------------------------------------------------
!     finished so return
!-----------------------------------------------------------------------
      return
      end
```
<p>and parameter file.</p>
```1.              ! domain height
1.              ! domain width
300.        ! Reynolds
10.             ! Peclet
0.              ! Injection start
0.              ! Injection end
10          ! cells in X
10          ! cells in Y
5000            ! max iterations
1.e-6       ! convergence criteria
0.6             ! P underrelaxation
0.6     ! U underrelaxation
0.6             ! V underrelaxation
1.              ! Scalar underrelaxation
```
<p>I believe by modifying the boundary conditions in UMomentum,Vmomentum, Pressure Correction and Scalar I can convert it to a lid driven cavity. as south is a known wall in the channel flow I based my edits on this fact when converting the east and west inlets to walls as such I used these as the correct boundary conditions.</p>
<p>for U Momentum function</p>
```! set boundary conditions for coefficients in the equations
! south
    aw(:,1)=0.
    ae(:,1)=0.
    as(:,1)=0.
    an(:,1)=0.
    ap(:,1)=1.
    su(:,1)=0.
! west
    aw(1,:)=0.
    ae(1,:)=0.
    as(1,:)=0.
    an(1,:)=0.
    ap(1,:)=1.
    su(1,:)=0.
! north
    aw(:,iNyUNodes)=0.
    ae(:,iNyUNodes)=0.
    as(:,iNyUNodes)=0.
    an(:,iNyUNodes)=0.
    ap(:,iNyUNodes)=1.
    su(:,iNyUNodes)=1.
! east
    aw(iNxUNodes,:)=0.
    ae(iNxUNodes,:)=0.
    as(iNxUNodes,:)=0.
    an(iNxUNodes,:)=0.
    ap(iNxUNodes,:)=1.
    su(iNxUNodes,:)=0.
```
<p>as there is no momentum in the Y direction I believe these to be the correct Boundary Conditions for VMomentum Function.</p>
```! south
    aw(:,1)=0.
    ae(:,1)=0.
    as(:,1)=0.
    an(:,1)=0.
    ap(:,1)=1.
    su(:,1)=0.
! west
    aw(1,:)=0.
    ae(1,:)=0.
    as(1,:)=0.
    an(1,:)=0.
    ap(1,:)=1.
    su(1,:)=0.
! north
    aw(:,iNyVNodes)=0.
    ae(:,iNyVNodes)=0.
    as(:,iNyVNodes)=0.
    an(:,iNyVNodes)=0.
    ap(:,iNyVNodes)=1.
    su(:,iNyVNodes)=0.
! east
    aw(iNxVNodes,:)=0.
    ae(iNxVNodes,:)=0.
    as(iNxVNodes,:)=0.
    an(iNxVNodes,:)=0.
    ap(iNxVNodes,:)=1.
    su(iNxVNodes,:)=0.
```

I am unsure how to manipulate the Pressure Corrections and Scalar Function's Boundary Conditions.

and the output of my code is either incorrect values or NaN with different attempts of the boundary conditions

Apologies for this being long, I wanted to be as thorough as possible as I cannot understand what is incorrect.

(it is known that the Fortran solver is validated as accurate, and my MATLAB solver) however these results do not matchup after my boundary edits.

Thanks in advance


问:
snappyHexMesh not meshing correct region
I created a simple model of cylindrical container to perform heat transfer simulation using FreeCAD for CAD modelling. I exported the model in stl format along with domain boundaries and faces. I have attached all the files here.

I have been trying to mesh the container domain. I am using cyclic blockMesh first to create a hex cylindrical block. Next I am using surfaceFeatures and snappyHexMesh to mesh the region inside the cylindrical container. To generate mesh in correct region I located a point using paraview that was clearly within the region where I want the mesh.

However, when I check the mesh in paraview, the mesh is generated for wrong region.
I tried different point locations but the result turns out to be the same.
I have attached the outputs of snappy and blockmesh as well.

Could someone help me understand why am I having this problem and how to solve it. Is it a bug in snappyHexMesh?

Thank you.

Here is the link to the files: Container_Simulation


问:
Conceptual doubt regarding 2D conjugate heat transfer modelling (COMSOL and Mathemtica)
I have been dealing with some conceptual flaws in my understanding of modelling, which I will elaborate herein. I am modelling conjugate heat transfer of a reciprocating fluid, which flows with velocity $u_{0} \sin(2 \pi f t)$, over a thick heated (heat input $Q$) solid of length $L$ and width $W$. I am initially modelling the problem in $2D$ in Mathematica, where I solve the coupled flow and energy equation (in both the domains) using the in-built fem approach. The non-dimensional energy equation I solve, looks like this:

\begin{equation}
\gamma Pe\left(\frac{\partial T}{\partial t}+(\vec{V}\cdot\nabla)T \right)=\nabla\cdot\left(\gamma\nabla T\right)
 \end{equation}

where $Re=u_0d/\nu$ is the Reinolds number, $Pe$ is the variable in space Peclet number

\begin{equation}
  Pe=\begin{cases}
   u_0d/\alpha_s, & \{x,y\}\in solid \\
   u_0d/\alpha_f, & \{x,y\}\in fluid \\
  \end{cases}
\end{equation}

Coefficient $\gamma$:

\begin{equation}
  \gamma=\begin{cases}
   1, & \{x,y\}\in solid \\
   k_f/k_s, & \{x,y\}\in fluid \\
  \end{cases}
\end{equation}
Since, its a $2D$ approach, my model only has the length $L \space (x-direction)$ and depth $e \space (y-direction)$ features. The width $W \space (z-direction)$ is not explicitly modelled. Hence, the heat input is supplied as a flux condition at the base:

$$q=\frac{Q}{L W}$$

Now, I model the same problem in the COMSOL simulation environment, again using a $2D$ approach. However, COMSOL asks me for a depth $d_z$ while setting up the physics, which I guess is the $W$ of my geometry. COMSOL solves the following equation:

As you can see, that the governing equation in COMSOL is multiplied by $d_z$ all over. The same is the case for its constant heat flux b.c., which looks like this

Hence, I supplied the same $q$ that I supplied to Mathematica model, while entering $W$ as the value of $d_z$ (its not 1m as in the picture, it should be 0.035m).

However, I cannot get the temperature evolution to match at all, there is a huge difference between Mathematica and COMSOL results even after grid and time-independence tests. The following plot shows the temperature history of a point in solid for the COMSOL and Mathematica model (all flow, material and geometry remain equivalent)

My questions are:

Is the problem in the amount of heat flux I am supplying to any one of the models ?

Let us assume from experiments, I know the $Q$ (in W) heat, which is supplied. Now imagine, I am modelling this system in 2D (in Mathematica), where my model only has the length ($L$) feature and not the width ($W$). In this scenario, I have to emulate the heat input as a flux ($q$) condition (Neumann value). How should I calculate this flux, to replicate the experiment ? Should it be $q = \frac{experimetal \space Q}{L∗1}$ or $q = \frac{experimental \space Q}{L∗W}$? Similar argument also follows for deciding the uavg value for the simulation as volumetric flow rate is known from experiments. Should it be $\frac{experimental \space flowate}{d*1}$ or $\frac{experimental \space flowate}{d*W}$ for the simulation.


问:
Delaunay-based isosurface extraction vs marching cubes
I recently tried the isosurface extraction algorithm provided by the C++ library CGAL. This is new to me. It is based on Delaunay triangulations.

I have some experience with the marching cubes, I tried to compare these two methods but it  is not easy to conclude. Regarding the speed, I have the impression that there's no winner, each method can beat the other, depending on the case. For example, the marching cubes algorithm is undoubtly faster for the trihyperboloid:

Regarding the quality of the mesh, I didn't observe some difference. I'm wondering whether there exist some examples for which the Delaunay-based method is better than the marching cubes. Would you know such examples? The Delaunay-based method has three input parameters (three bounds, an angle, a distance and a radius), while the marching cubes only requires the choice of a grid. So I believe it's possible that the Delaunay-based method can achieve something better in some cases.


问:
Solving 2D Heat Equation with Input using central finite difference method
So I have implemented central finite difference method for solving the 2D heat equation. When I leave all initial and boundary conditions as 0s, but apply an input uniform across the entire space or grid, I am only noticing temperature response at the boundary. I am not sure exactly what is causing this.

The code is here: https://colab.research.google.com/drive/1YCh_DDrw0uZPLvrWKU2hcw6lSe7cQiRS

But for a quick example of what this looks like:

Would really appreciate any help with looking at the code to see what's causing this. Is there any chance that my algorithm is unstable and I need to increase the discretization points or something?