How to remove rotational modes (i.e., Null space) in sphere?

[gthyagi]

Here are the details of null space removal from the benchmark paper.

https://gmd.copernicus.org/articles/14/1899/2021/#:~:text=Solving%20the%20linear%20system%20and%20dealing%20with%20null%20spaces

Do we also need to remove zero modes from the approximate solution at every iteration?

[knepley]

Yes. If you set the nullspace of the Mat object, the solver will do it automatically.

[lmoresi]

Your previous suggestion was MG with SVD as the coarse solve, what's your view on that @knepley ?

[knepley]

If you can use MG, I think that is the best way to go. You might still need it for the outer Krylov method, depending on exactly what smoother/prolongator choice you make, but it will pick it up if it is attached to the fine matrix.

[gthyagi]

I am using following Stokes settings

# Stokes settings

stokes.tolerance = stokes_tol
stokes.petsc_options["ksp_monitor"] = None

stokes.petsc_options["snes_type"] = "newtonls"
stokes.petsc_options["ksp_type"] = "fgmres"

# stokes.petsc_options.setValue("fieldsplit_velocity_pc_type", "mg")
stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_type", "kaskade")
stokes.petsc_options.setValue("fieldsplit_velocity_pc_mg_cycle_type", "w")

stokes.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd"

stokes.petsc_options[f"fieldsplit_velocity_ksp_type"] = "fcg"
stokes.petsc_options[f"fieldsplit_velocity_mg_levels_ksp_type"] = "chebyshev"
stokes.petsc_options[f"fieldsplit_velocity_mg_levels_ksp_max_it"] = 5
stokes.petsc_options[f"fieldsplit_velocity_mg_levels_ksp_converged_maxits"] = None

# # gasm is super-fast ... but mg seems to be bulletproof
# # gamg is toughest wrt viscosity
# stokes.petsc_options.setValue("fieldsplit_pressure_pc_type", "gamg")
# stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_type", "additive")
# stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_cycle_type", "v")

# mg, multiplicative - very robust ... similar to gamg, additive
stokes.petsc_options.setValue("fieldsplit_pressure_pc_type", "mg")
stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_type", "multiplicative")
stokes.petsc_options.setValue("fieldsplit_pressure_pc_mg_cycle_type", "v")

I am subtracting following component to removed the nullspace in velocity.
Annulus:

# Null space in velocity (constant v_theta) expressed in x,y coordinates
v_theta_fn_xy = r_uw * mesh.CoordinateSystem.rRotN.T * sympy.Matrix((0,1))
print(v_theta_fn_xy)

$$ \begin{bmatrix} -y \\ x \end{bmatrix} $$

# Null space evaluation
I0 = uw.maths.Integral(mesh, v_theta_fn_xy.dot(v_uw.sym))
norm = I0.evaluate()
I0.fn = v_theta_fn_xy.dot(v_theta_fn_xy)
vnorm = I0.evaluate()

# print(norm/vnorm, vnorm)

with mesh.access(v_uw):
    dv = uw.function.evaluate(norm * v_theta_fn_xy, v_uw.coords) / vnorm
    v_uw.data[...] -= dv 

Spherical Shell:

# Null space in velocity expressed in x,y,z coordinates
v_theta_phi_fn_xyz = sympy.Matrix(((0,1,1), (-1,0,1), (-1,-1,0))) * mesh.CoordinateSystem.N.T
print(v_theta_phi_fn_xyz)

$$ \begin{bmatrix} y+z \\ -x+z \\ -x-y \end{bmatrix} $$

# Null space evaluation

I0 = uw.maths.Integral(mesh, v_theta_phi_fn_xyz.dot(v_uw.sym))
norm = I0.evaluate()

I0.fn = v_theta_phi_fn_xyz.dot(v_theta_phi_fn_xyz)
vnorm = I0.evaluate()
# print(norm/vnorm, vnorm)

with mesh.access(v_uw):
    dv = uw.function.evaluate(norm * v_theta_phi_fn_xyz, v_uw.coords) / vnorm
    v_uw.data[...] -= dv 

I have taken the rotation matrix for spherical sphell from here:
chrome-extension://efaidnbmnnnibpcajpcglclefindmkaj/https://tigerprints.clemson.edu/cgi/viewcontent.cgi?article=5105&context=all_theses

Is the nullspace removal for the spherical shell correct?
Nullspace removal in Annulus case is implemented/suggested by @lmoresi and I am extending it to spherical shell.

[lmoresi]

In the spherical examples, I removed three orthogonal rotation terms sequentially and didn't worry about any mathematical fanciness. I think this is the same as your approach above.

orientation_wrt_z = sympy.atan2(y + 1.0e-10, x + 1.0e-10)
v_rbm_z_x = -ra * sympy.sin(orientation_wrt_z)
v_rbm_z_y = ra * sympy.cos(orientation_wrt_z)
v_rbm_z = sympy.Matrix([v_rbm_z_x, v_rbm_z_y, 0]).T

orientation_wrt_x = sympy.atan2(z + 1.0e-10, y + 1.0e-10)
v_rbm_x_y = -ra * sympy.sin(orientation_wrt_x)
v_rbm_x_z = ra * sympy.cos(orientation_wrt_x)
v_rbm_x = sympy.Matrix([0, v_rbm_x_y, v_rbm_x_z]).T

orientation_wrt_y = sympy.atan2(z + 1.0e-10, x + 1.0e-10)
v_rbm_y_x = -ra * sympy.sin(orientation_wrt_y)
v_rbm_y_z = ra * sympy.cos(orientation_wrt_y)
v_rbm_y = sympy.Matrix([v_rbm_y_x, 0, v_rbm_y_z]).T