Stokes3D.py Incurs High Error after Changing the Coefficients of the Periodic Flow on the X-Y Plane

[MrMonotreme]

When one changes uex = sin(2*y)*(1-z**2) to uex = sin(2.05*y)*(1-z**2) in the manufactured solution, the error for the x-component of the flow velocity increases from 3.9229e-14 to 2.672e0.

Would I need to change some of the parameters when creating the Fourier basis for the periodic flow on the x-axis?

Furthermore, would it be possible to perform a singular value decomposition on the non-exact solution produced by Shenfun to determine which basis functions preserve the greatest variance of the fluid dynamics. Would the following code suffice?

r = 5
for dimension in range(0, 3):
            u_decomp, s_decomp, vh_decomp = np.linalg.svd(u_[dimension, :, :, :])
            print(u_decomp.shape, s_decomp.shape, vh_decomp.shape)
            for i in range(n):
                plt.figure()
                plt.yscale('log')
                plt.scatter(np.arange(n), s_decomp[i])
                plt.savefig("./sigma/" + str(dimension) + "_" + str(i) + "_variance_decomp.png")
            plt.figure()
            plt.yscale('linear')
            for i in range(0, r):
                plt.plot(np.arange(n), u_decomp[i])
            plt.savefig(str(dimension) + "_variance_modes.png")

variance_decomposition.png

variance_modes.png
The resulting modes are quite jumbled as I am having trouble separating the various axes of the block matrix M.

meshgrid.png
Here is my pictorial understanding of the local mesh X. However, I don't have an intuitive grasp on block matrix M.

[mikaem]

Would I need to change some of the parameters when creating the Fourier basis for the periodic flow on the x-axis?

The manufactured solutions should be periodic, otherwise a large error is incurred. If you absolutely need to use sin(2.05*y), then the domain must be changed accordingly.

Furthermore, would it be possible to perform a singular value decomposition on the non-exact solution produced by Shenfun to determine which basis functions preserve the greatest variance of the fluid dynamics. Would the following code suffice?

Not sure what your true objective is here? If you are looking for which basis functions that have the greatest contributions to the solution I think you should look at the coefficients up_hat and not the solution in real space up.

Regarding the block matrix M, remember that this matrix represents the Legendre/Chebyshev part of the problem, and there is one matrix M for each Fourier coefficient. I have tried to explain this in the Stokes demo

[MrMonotreme]

The manufactured solutions should be periodic, otherwise a large error is incurred. If you absolutely need to use sin(2.05*y), then the domain must be changed accordingly.

Would it be possible to use uex= sin(2.05*y)*(1-z**2) and uey= sin(2*x)*(1-z**2)? I characterized the periodic solutions on the x-y plane over a domain from 0 to 2π and set N, the number of quadrature points, to 400.
K0 = Basis(400, 'Fourier', dtype='D', domain=(0, 2*np.pi))
I received an error of 6.7283e+01 on the x axis, which was higher than 4.1541, the error incurred when I set N to 15. Are domains of fewer quadrature points preferable to those of more quadrature points when characterizing periodic flows that are minutely out of phase?

Not sure what your true objective is here? If you are looking for which basis functions that have the greatest contributions to the solution I think you should look at the coefficients up_hat and not the solution in real space up.

I'm trying to determine which preserve the greatest variance of the periodic flows in the x-y plane. Ultimately, I hope to reduce the number of modes used to approximate the dynamics of the Navier-Stokes equations. Does this functionality already exist in shenfun? I may have looked in the wrong place of the documentation.

Is shenfun able to characterize the behavior of fluids, which are prone to cavitations?

Thank you for your support. This library has saved me countless hours.

[mikaem]

Would it be possible to use uex= sin(2.05y)(1-z2) and uey= sin(2x)(1-z2)?

Not without changing the domain. You would need to use for example

K1 = Basis(N[1], 'Fourier', dtype='d', domain=(0, 2*np.pi/2.05))

or perhaps domain=(0, 4*np.pi/2.05). This is just because sin(2.05*y) then will be periodic on the domain. The manufactured solution must be periodic on the domain since the Fourier basis functions are periodic on the domain. A non periodic solution will always have a poor approximation with Fourier bases.

Reducing the number of modes is easy, there are already functions available for this. See refine and assign.

About the cavitations, could you please be more specific?

[MrMonotreme]

I'm sorry for the late response. I have been able to successfully form a reduced order model of the fluid dynamics using domain=(0, 4*np.pi/2.05) and refine the model to fewer modes.

As a separate project, I've been trying to apply this to the mhd solver in your gist.

For tgmhd.py on line 33, I receive

AttributeError: 'TensorProductSpace' object has no attribute 'local_shape'

K2 = np.zeros(T.local_shape(True), dtype=float) where T = TensorProductSpace(comm, (V0, V1, V2), **{'planner_effort': 'FFTW_MEASURE'})
Is local_shape still available for TensorProductSpaces? I checked the history for tensorproductspace.py and the documentation.

Should I be using shape or global_shape() instead?

[mikaem]

It should be just shape. Makes no difference if you use only one cpu.

[MrMonotreme]

I was able to run the program without a crash by using the following.

K2 = np.zeros(T.shape(), dtype=float)
for i in range(dim):
    K2 += K[i]*K[i]

K_over_K2 = np.zeros((3, 32, 32, 17), dtype=float) # TV = vector
for i in range(dim):
    K_over_K2[i] = K[i] / np.where(K2 == 0, 1, K2)

Would this be sufficient for fixing the error?
I received the following output.

K 0.12456540817722132 ERR 2.2132295995902496e-13 B 0.12463776214378505 ERR 7.85052578500256e-13
K 0.12456540817673264 ERR -2.673555821175455e-13 B 0.12463776214353375 ERR 5.337535968763518e-13

When calculating the pressure gradient dU[i] -= P_hat*K[i], could I easily form a new variable which calculates the Laplacian of the pressure?
Would something as simple as multiplying by the local wavenumbers twice P_hat*K[i]*K[i] calculate the Laplacian of the pressure?

[MrMonotreme]

When running fh_hat.refine((5, 5, 5)) repeatedly, I occasionally receive a traceback of

Traceback (most recent call last):
  File "/home/username/anaconda3/envs/neppu/lib/python3.8/site-packages/shenfun/forms/arguments.py", line 1129, in refine
    base = space.bases[axes[0]]
  File "/home/username/anaconda3/envs/neppu/lib/python3.8/site-packages/shenfun/tensorproductspace.py", line 906, in __getattr__
    assert name not in ('bases',)

Does the refine method occasionally fail for Functions of dimensionality greater than one? The intended functionality f_hat.refine((5, 5, 5)) works when called before up_hat= M.solve(fh_hat, constraints=((3, 0, 0), (3, N[2]-1, 0))).

Is there a general method for performing something similar to https://youtu.be/X5GhhjpX0ao?t=2145? Or, could someone form a new model of the fluid dynamics after solving the original system via the blockmatrix? Could the modes of this new model be chosen to preserve the greatest variance of the original system using a technique such as singular value decomposition?