Fixed mssing operators in viscous force

repepo · May 24, 2024 · 62426d7 · 62426d7
1 parent e6fe4fe
commit 62426d7
Show file tree

Hide file tree

Showing 2 changed files with 24 additions and 26 deletions.
diff --git a/bin/operators_new.py b/bin/operators_new.py
@@ -153,20 +153,20 @@ def viscous_diffusion(l, section, component, offdiag):  # ----------------------
 
                 if par.variable_viscosity:
 
-                    out = L * ( -L*(l+2)*(l-1)*r0_vsc0_D0_u - (L+2)*r1_vsc0_lho1_D0_u - 2 * (L+1) * r1_vsc1_D0_u
-                                -2*(L-1)*(r2_vsc0_lho2_D0_u + r2_vsc1_lho1_D0_u) - (l+2)*(l-1)*r2_vsc2_D0_u
-                                + r3_vsc0_lho3_D0_u + 2*r3_vsc1_lho2_D0_u + r3_vsc2_lho1_D0_u
+                    out = L * ( L*(2-L)* r0_vsc0_D0_u - (L+2)* r1_vsc0_lho1_D0_u - 2*(L+1)* r1_vsc1_D0_u
+                                -2*(L-1)*( r2_vsc0_lho2_D0_u + r2_vsc1_lho1_D0_u ) + (2-L)* r2_vsc2_D0_u
+                                + r3_vsc0_lho3_D0_u + 2* r3_vsc1_lho2_D0_u + r3_vsc2_lho1_D0_u
 
-                                + (2-L)*r2_vsc0_lho1_D1_u + 6 * (r3_vsc0_lho2_D1_u + r3_vsc1_lho1_D1_u)
-                                + r4_vsc0_lho3_D1_u + 2*r4_vsc1_lho2_D1_u + r4_vsc2_lho1_D1_u
+                                + (2-L)* r2_vsc0_lho1_D1_u + 6* ( r3_vsc0_lho2_D1_u + r3_vsc1_lho1_D1_u )
+                                + r4_vsc0_lho3_D1_u + 2* r4_vsc1_lho2_D1_u + r4_vsc2_lho1_D1_u + 2*(L+1)* r2_vsc1_D1_u
 
-                                + 2*L*r2_vsc0_D2_u + 5*r3_vsc0_lho1_D2_u - 4*r3_vsc1_D2_u + 2*r4_vsc0_lho2_D2_u
-                                + 2*r4_vsc1_lho1_D2_u - r4_vsc2_D2_u
+                                + 2*L* r2_vsc0_D2_u + 5* r3_vsc0_lho1_D2_u - 4* r3_vsc1_D2_u + 2* r4_vsc0_lho2_D2_u
+                                + 2* r4_vsc1_lho1_D2_u - r4_vsc2_D2_u
 
-                                -4*r3_vsc0_D3_u + r4_vsc0_lho1_D3_u - 2*r4_vsc1_D3_u
+                                -4* r3_vsc0_D3_u + r4_vsc0_lho1_D3_u - 2* r4_vsc1_D3_u
+
+                                - r4_vsc0_D4_u )
 
-                                - r4_vsc0_D4_u
-                                )
             else:
 
                 if (par.magnetic == 1 and par.B0 == 'dipole'):
@@ -182,13 +182,10 @@ def viscous_diffusion(l, section, component, offdiag):  # ----------------------
                             +r2_D2_v)
 
                 if par.variable_viscosity:
-
-                    out = L * ( -L*r0_vsc0_D0_v - 3*r1_vsc0_lho1_D0_v - r2_vsc1_lho1_D0_v
-
-                            + 2*r1_vsc0_D1_v - r2_vsc0_lho1_D1_v + r2_vsc1_D1_v
-
-                            + r2_vsc0_D2_v
-                            )
+                    out = L * ( -L* r0_vsc0_D0_v - 3* r1_vsc0_lho1_D0_v - r2_vsc1_lho1_D0_v - r2_vsc0_lho2_D0_v - r1_vsc1_D0_v
+                            + 2* r1_vsc0_D1_v - r2_vsc0_lho1_D1_v + r2_vsc1_D1_v
+                            + r2_vsc0_D2_v )
+
             else:
                 if (par.magnetic == 1 and par.B0 == 'dipole'):
                     out = L*( -L*r3_D0_v + 2*r4_D1_v + r5_D2_v )                          # r5* r.1curl( nabla^2 u )

diff --git a/bin/submatrices_new.py b/bin/submatrices_new.py
@@ -200,13 +200,13 @@ def main(ncpus):
 
             # And even more viscous diffusion,
             if par.variable_viscosity:
-                labl_u += ['r0_vsc0_D0', 'r1_vsc0_lho1_D0', 'r1_vsc1_D0', 'r2_vsc0_D2',
-                        'r2_vsc0_lho1_D1', 'r2_vsc0_lho2_D0', 'r2_vsc1_lho1_D0',
-                        'r2_vsc2_D0', 'r3_vsc0_D3', 'r3_vsc0_lho1_D2', 'r3_vsc0_lho2_D1',
-                        'r3_vsc0_lho3_D0', 'r3_vsc1_D2', 'r3_vsc1_lho1_D1', 'r3_vsc1_lho2_D0',
-                        'r3_vsc2_lho1_D0', 'r4_vsc0_D4', 'r4_vsc0_lho1_D3', 'r4_vsc0_lho2_D2',
-                        'r4_vsc0_lho3_D1', 'r4_vsc1_D3', 'r4_vsc1_lho1_D2', 'r4_vsc1_lho2_D1',
-                        'r4_vsc2_D2', 'r4_vsc2_lho1_D1']
+                labl_u += ['r0_vsc0_D0'     , 'r1_vsc0_lho1_D0', 'r1_vsc1_D0'     , 'r2_vsc0_D2'     ,
+                           'r2_vsc0_lho1_D1', 'r2_vsc0_lho2_D0', 'r2_vsc1_lho1_D0', 'r2_vsc1_D1'     ,
+                           'r2_vsc2_D0'     , 'r3_vsc0_D3'     , 'r3_vsc0_lho1_D2', 'r3_vsc0_lho2_D1',
+                           'r3_vsc0_lho3_D0', 'r3_vsc1_D2'     , 'r3_vsc1_lho1_D1', 'r3_vsc1_lho2_D0',
+                           'r3_vsc2_lho1_D0', 'r4_vsc0_D4'     , 'r4_vsc0_lho1_D3', 'r4_vsc0_lho2_D2',
+                           'r4_vsc0_lho3_D1', 'r4_vsc1_D3'     , 'r4_vsc1_lho1_D2', 'r4_vsc1_lho2_D1',
+                           'r4_vsc2_D2'     , 'r4_vsc2_lho1_D1' ]
 
         if par.magnetic == 1 :
             # add Lorentz force
@@ -253,8 +253,9 @@ def main(ncpus):
             labl_v += [ 'r1_lho1_D0', 'r2_lho2_D0', 'r2_lho1_D1' ]
 
             if par.variable_viscosity:
-                labl_v += ['r0_vsc0_D0', 'r1_vsc0_D1', 'r1_vsc0_lho1_D0', 'r2_vsc0_D2',
-                           'r2_vsc0_lho1_D1', 'r2_vsc1_D1', 'r2_vsc1_lho1_D0']
+
+                labl_v += ['r0_vsc0_D0'     , 'r1_vsc0_D1', 'r1_vsc0_lho1_D0', 'r2_vsc0_D2'     ,
+                           'r2_vsc0_lho1_D1', 'r2_vsc1_D1', 'r2_vsc1_lho1_D0', 'r2_vsc0_lho2_D0', 'r1_vsc1_D0']
 
         if par.magnetic == 1 :
             # Lorentz force