Try to rescue periodic ODE solver in latest changes by distinguishing…

… between periodic BC modifications required for matrices on LHS or RHS in LHS.x = RHS. The ODE solver in the script sets RHS[end] = 0.0 so second version of mass matrix is avoided. However, L_matrix_with_bc appears to be singular in some examples, suggesting a continuing problem.
mabarnes · Sep 12, 2024 · b76c05d · b76c05d
1 parent 734ab81
commit b76c05d
Showing 1 changed file with 8 additions and 8 deletions.
diff --git a/moment_kinetics/src/gauss_legendre.jl b/moment_kinetics/src/gauss_legendre.jl
@@ -122,13 +122,13 @@ function setup_gausslegendre_pseudospectral(coord; collision_operator_dim=true)
 
     dirichlet_bc = (coord.bc in ["zero", "constant"]) # and further options in future
     periodic_bc = (coord.bc == "periodic")
-    setup_global_weak_form_matrix!(mass_matrix, lobatto, radau, coord, "M"; periodic_bc=periodic_bc)
-    setup_global_weak_form_matrix!(K_matrix, lobatto, radau, coord, "K_with_BC_terms"; periodic_bc=periodic_bc)
+    setup_global_weak_form_matrix!(mass_matrix, lobatto, radau, coord, "M"; periodic_bc_lhs=periodic_bc)
+    setup_global_weak_form_matrix!(K_matrix, lobatto, radau, coord, "K_with_BC_terms"; periodic_bc_rhs=periodic_bc)
     setup_global_weak_form_matrix!(L_matrix, lobatto, radau, coord, "L_with_BC_terms")
     mass_matrix_lu = lu(sparse(mass_matrix))
-    if dirichlet_bc #|| periodic_bc
+    if dirichlet_bc || periodic_bc
         L_matrix_with_bc = allocate_float(coord.n,coord.n)
-        setup_global_weak_form_matrix!(L_matrix_with_bc, lobatto, radau, coord, "L", dirichlet_bc=dirichlet_bc, periodic_bc=periodic_bc)
+        setup_global_weak_form_matrix!(L_matrix_with_bc, lobatto, radau, coord, "L", dirichlet_bc=dirichlet_bc, periodic_bc_lhs=periodic_bc)
         L_matrix_with_bc = sparse(L_matrix_with_bc )
         L_matrix_lu = lu(sparse(L_matrix_with_bc))
     else
@@ -880,7 +880,7 @@ is supported with boundary conditions.
 function setup_global_weak_form_matrix!(QQ_global::Array{mk_float,2},
                                lobatto::gausslegendre_base_info,
                                radau::gausslegendre_base_info, 
-                               coord,option; dirichlet_bc=false, periodic_bc=false)
+                               coord,option; dirichlet_bc=false, periodic_bc_lhs=false, periodic_bc_rhs=false)
     QQ_j = allocate_float(coord.ngrid,coord.ngrid)
 
     ngrid = coord.ngrid
@@ -915,19 +915,19 @@ function setup_global_weak_form_matrix!(QQ_global::Array{mk_float,2},
         end
         # requires RHS vector b[1],b[end] = boundary values
     end
-    if periodic_bc
+    if periodic_bc_lhs || periodic_bc_rhs
         # Make periodic boundary condition by modifying elements of matrix for duplicate point
         # add assembly contribution to lower endpoint from upper endpoint
         j = coord.nelement_local
         get_QQ_local!(QQ_j,j,lobatto,radau,coord,option)
-        iminl = imin[j] - (coord.nelement_local > 1)
+        iminl = imin[j] - mk_int(coord.nelement_local > 1)
         imaxl = imax[j]
         QQ_global[1,iminl:imaxl] .+= QQ_j[end,:]
         # Enforce continuity at the periodic boundary
         # All-zero row in RHS matrix sets last element of `b` vector in
         # `mass_matrix.x = b` to zero
         QQ_global[end,:] .= 0.0
-        if option == "M"
+        if periodic_bc_lhs
             # enforce periodicity `x[1] = x[end]` using the last row of the
             # `mass_matrix.x = b` system.
             QQ_global[end,1] = 1.0