Test Stokes Identities

[alexfikl]

Add tests for some identities satisfied by the Stokeslet and the Stresslet on arbitrary geometries. They're given on page 2 of
https://services.math.duke.edu/~leili/teaching/duke/math660s18/lectures/lec25.pdf
or in Pozrikidis' tome of knowledge if available.

[alexfikl]

Convergence results in 2D on a 5-armed starfish for the Stokeslet

h	5rd order	EOC	3rd order	EOC
0.0625	7.347227267334338e-12		1.155847485638326e-09
0.03125	6.995540192862465e-12	00.07	9.154592966619713e-10	0.33
0.015625	4.662484448865704e-15	10.55	3.639546560164975e-12	7.97
0.01041666	2.742219800926962e-17	12.66	1.136946078694970e-13	8.54
0.0078125	6.739665360935961e-19	12.88	7.701738440580504e-15	9.35

and the Stresslet

h	5th order	EOC	3rd order	EOC
0.0625	7.841613331065972e-08		5.169261460949249e-06
0.03125	7.335107515375984e-08	00.09	3.319709609159663e-06	0.63
0.015625	1.169692803219092e-10	09.29	4.870965283109148e-08	6.09
0.01041666	1.323785820170431e-12	11.05	2.415487781430483e-09	7.40
0.0078125	4.911512430291654e-14	11.45	2.452037085419880e-10	7.95

No idea where the orders are coming from here. Looks like a quadrature 2n + 1 order (?), but I don't see why that error would dominate.

EDIT: Nevermind, managed to make a stupid mistake in the norm. Used norm(x.dot(x)) instead of norm(sqrt(x.dot(x))), so got an extra square in there. The orders are 3-4 and 5-6 now, as expected.

[alexfikl]

Convergence in 3D on a sphere with target_order = qbx_order = 3 for the Stokeslet

h	4x	EOC	8x	EOC
1.2846885921500641	0.0001710970985557689		3.148996157356721e-07
0.6597867118634438	6.563532962691179e-05	1.43	2.864323608373111e-09	7.05
0.3307289701916807	2.716344909532715e-05	1.27	5.346893474359467e-10	2.43

and the Stresslet

h	4x	EOC	8x	EOC
1.2846885921500641	0.0673876952443450		7.04018124819349e-05
0.6597867118634438	0.0689589261268387	-0.03	1.38161221816446e-06	5.89
0.3307289701916807	0.0651196246940797	+0.08	4.48320362983306e-07	1.62

and for target_order = qbx_order = 5 for the Stresslet (since at 3rd order it seemed to not wor at all)

h	4x	EOC
1.3015230810921776	0.005921827361342184
0.6593342656530706	0.004444005389353953	0.42
0.33051732288684144	0.004126663149293894	0.10

At 8x it seems to converge for a bit, but it bottoms out fairly quickly. At 4x its at best first order.

[alexfikl]

Just because it's so pretty, this is the error in the Stokeslet for target_order = qbx_order = 3 and 4x oversampling and the finest refinement.

[inducer]

Thanks for doing this! So overall this seems to match results from #45 pretty well (IIUC)... something is causing accuracy to degrade severely in 3D for Stokes. In a way, I'm glad it's not a solve-only phenomenon. Building on the conjecture that the geometry derivatives are inaccurate: the QBX radius and center placement are derived from geometry derivatives, and if these are off, then I could very easily see things coming off the rails. Do things improve if you replace those with analytic/constant quantities?

(Edit: It's still puzzling though that things like the Laplace double layer wouldn't be affected... but the Stokes kernels are "more directional"? Don't know.)

[alexfikl]

Building on the conjecture that the geometry derivatives are inaccurate: the QBX radius and center placement are derived from geometry derivatives, and if these are off, then I could very easily see things coming off the rails. Do things improve if you replace those with analytic/constant quantities?

Ah, completely forgot about those! I can easily see them actually having a bigger effect, as opposed to the area elements or other things in there. I'll try that out.

[alexfikl]

Ok, so time for another data dump with recent experiments. All of these are for the Stokeslet identity test with target_order = qbx_order = 3 and 4x oversampling. Some of the numbers are slightly different than before because it's using pytential.norm to compute the relative errors.

First off, using analytical normals + a constant expansion radius (basically h_max * Ones()) for all elements seems to help a bit.

h	Baseline	EOC	Normals	EOC
1.28468859e+00	4.90769308e-04		2.60256347e-03
6.59786712e-01	1.23417237e-04	2.07	2.66794923e-04	3.42
3.30728970e-01	3.69521709e-05	1.75	3.71973115e-05	2.85
1.65394094e-01	1.54490224e-05	1.26	1.28638136e-05	1.53
8.26980033e-02	7.17999555e-06	1.11	5.59909245e-06	1.20

Next, I tried playing with the fudge factor for the expansion radii. By default it was set to 1 in 3D, so I tried bumping it up to 1.25 and 1.5. This should still give h-convergence, right?

h	1.25	EOC	1.5	EOC	1.25 + Normals	EOC
1.28468859e+00	1.42968485e-04		5.87162861e-05		2.25151878e-03
6.59786712e-01	2.47939804e-05	2.63	7.55798051e-06	3.08	1.95342839e-04	3.67
3.30728970e-01	5.37776584e-06	2.21	9.41696442e-07	3.02	1.23640652e-05	4.00
1.65394094e-01	1.96232164e-06	1.45	2.92240655e-07	1.69	1.79200673e-06	2.79
8.26980033e-02	8.54446355e-07	1.20	1.18752749e-07	1.30	6.53734352e-07	1.45

Also tried some of these combinations on the Stresslet identity with the same parameters. They seem to bottom out to first order a lot faster, but definitely better than the zero-order from #108 (comment).

h	1.25 + Normals	EOC
1.28468859e+00	1.23881229e-02
6.59786712e-01	2.04313813e-03	2.70
3.30728970e-01	2.83060944e-04	2.86
1.65394094e-01	9.14596971e-05	1.63
8.26980033e-02	6.40306872e-05	0.51

[isuruf]

I tried with my branch with StokesletWrapperTornberg and fmm_order=10. Here is what I get for Stokeslet,

h              | Error          | Running EOC
---------------+----------------+-------------
1.30152308e+00 | 3.62742818e-02 |            
6.59334266e-01 | 1.66108442e-02 | 1.15       
3.30517323e-01 | 4.54252051e-03 | 1.88       
1.65361804e-01 | 8.80006450e-04 | 2.37       
8.26937755e-02 | 1.38277785e-04 | 2.67  

[alexfikl]

I tried with my branch with StokesletWrapperTornberg and fmm_order=10. Here is what I get for Stokeslet,

@isuruf What's the expected order for that? I seem to remember Andreas saying that has additional target derivatives, so it needs bumping the qbx order to get back to the same as the naive version (?).

[isuruf]

It has one target derivative and one source derivative. So, I would expect 2 orders reduced.

[isuruf]

It's based off of https://doi.org/10.1016/j.jcp.2007.06.029.

[alexfikl]

After our discussion today, added a few more data points here. First, using target_order = 5, qbx_order = 3 and 4x oversampling. On a sphere / spheroid and running

python test_stokes.py 'test_stokeslet_identity(_acf, partial(eid.SpheroidTestCase, resolutions=[0, 1, 2, 3]), visualize=True)'

h	Sphere	EOC	h	Spheroid	EOC
1.30152308e+00	5.28414578e-06		1.30152308e+00	1.13522864e-05
6.59334266e-01	2.14063344e-07	4.71	6.59334266e-01	2.18712820e-07	5.81
3.30517323e-01	3.98685216e-08	2.43	3.30517323e-01	2.27535794e-08	3.28
1.65361804e-01	1.34397959e-08	1.57	1.65361804e-01	8.70470303e-09	1.39

Probably important: using a fudge_factor for the expansion radii of 1.25 vs the default 0.5, so quite a lot larger.

[alexfikl]

@inducer This should be ready for a look now. It just adds two extra test cases for Stokes. Some takeaways:

• with the default dim_fudge_factor = 0.5, it requires 8x oversampling to look like it's converging
• bumping to dim_fudge_factor = 1.0 (like in 2D), it converges with 4x oversampling (most of the results above are in this regime).

[inducer]

Thanks! LGTM, just two minor nits, then it's good to go.

[inducer]

• Name arguments to designate one as reference?
• These look like two functions mashed into one with an if?

[alexfikl]

These look like two functions mashed into one with an if?

They were wearing a trenchcoat! How could you tell? 🤦

Cleaned up in 45efe7b.

[inducer]

Thanks!