To run the main Rosenbrock time stepper function (https://github.com/gridapapps/GridapGeosciences.jl/blob/master/src/ShallowWaterRosenbrock.jl#L8) from a parallel driver (while the parallel driver does not exist yet, we can take the Full Newton driver as a reference). This is likely to fail, as some of the lines in the Rosenbrock solver may not be prepared to be run from both environments (we now that it is prepared for serial environments, but didnt check in a parallel environment). Once we identify where it fails, we can try to understand why, and solve it, etc.

The good news is that we have the sequential driver to compare results against. These are our target results. It is not building something from scratch.

@amartinhuertas can provide support along the way

Note that I have partially done this work. For example, the rosenbrock time stepper no longer calls lu, lu!, ldiv! ; these are Julia interfaces. Instead, I replaced these calls by their counterparts in Gridap solver interfaces (That GridapPETSc implements). See, e.g., https://github.com/gridapapps/GridapGeosciences.jl/blob/master/src/ShallowWaterRosenbrock.jl#L121

This is working in the serial case, as, e.g., jacobian_matrix_solver, defaults to lu https://github.com/gridapapps/GridapGeosciences.jl/blob/master/src/ShallowWaterRosenbrock.jl#L78

I would start running the parallel program on just a single MPI task. Once we know it works, I would move to more than a single task.

see #35 (comment)

The failing tests arise with the following upgraded versions of GridapGeosiences dependencies

GridapDistributed v0.2.6
PartitionedArrays v0.2.10
GridapP4est 0.1.3
Gridap 0.17.13

After reading some comments in the code, it seems we have some issues to discuss:

1. The @law for dispersion is not working. We should at least try to identify which is the problem.
2. The doubt abut why a minus in the Neumann term. I think the answer is simple, since the problem is nonlinear, we are providing the solver the residual. The Neumann term is one of these terms, and in the residual it comes with a minus, right?

On the other hand, there are some additional developments:

1. Create tests with a method of manufactured solutions that will tell us whether the implementation is correct.
2. Implement the transient version of the driver using GridapODEs.

With regard to profiling:

1. Compare the computation time of the integration/assembly of Gridap vs FEMPAR. (Take into account that the first run involves JIT, consider a one warm-up solve first.)
2. Consider using a cheaper iterative e.g. multigrid solver or a PETSc solver. Right now, we are using the Julia \, not sure how efficient it is going to be.

What else?

Following the description here I ran JET on Julia 1..8.0-DEV with

JET.report_file("test/Williamson2ThetaMethodFullNewtonTests.jl")

which generated the attached report.

JET_Williamson2ThetaMethodFullNewtonTests.txt
.

To run @report_opt I had to modify the test from here to look like

module Williamsom2ThetaMethodFullNewtonTests

using Test
using Gridap
using GridapGeosciences
using SparseMatricesCSR

# Solves the steady state Williamson2 test case for the shallow water equations on a sphere
# of physical radius 6371220m. Involves a modified coriolis term that exactly balances
# the potential gradient term to achieve a steady state
# reference:
# D. L. Williamson, J. B. Drake, J. J.HackRüdiger Jakob, P. N.Swarztrauber, (1992)
# J Comp. Phys. 102 211-224

include("Williamson2InitialConditions.jl")

function test()
  l2_err_u = [0.011370921987771046 , 0.002991229234176333 ]
  l2_err_h = [0.005606685579166809, 0.001458107077200681 ]

  order=1
  degree=4
  θ=0.5
  for i in 1:2
    n      = 2*2^i
    nstep  = 5*n
    Uc     = sqrt(g*H₀)
    dx     = 2.0*π*rₑ/(4*n)
    dt     = 0.25*dx/Uc
    println("timestep: ", dt)   # gravity wave time step
    T      = dt*nstep
    τ      = dt/2
    model = CubedSphereDiscreteModel(n; radius=rₑ)
    hf, uf = shallow_water_theta_method_full_newton_time_stepper(model, order, degree,
                                                                 h₀, u₀, f₀, topography, g, θ, T, nstep, τ;
                                                                 write_solution=false,
                                                                 write_solution_freq=5,
                                                                 write_diagnostics=true,
                                                                 write_diagnostics_freq=1,
                                                                 dump_diagnostics_on_screen=true)

    Ω     = Triangulation(model)
    dΩ    = Measure(Ω, degree)
    hc    = CellField(h₀, Ω)
    e     = h₀-hf
    err_h = sqrt(sum(∫(e⋅e)*dΩ))/sqrt(sum(∫(hc⋅hc)*dΩ))
    uc    = CellField(u₀, Ω)
    e     = u₀-uf
    err_u = sqrt(sum(∫(e⋅e)*dΩ))/sqrt(sum(∫(uc⋅uc)*dΩ))
    println("n=", n, ",\terr_u: ", err_u, ",\terr_h: ", err_h)

    #@test abs(err_u - l2_err_u[i]) < 10.0^-12
    #@test abs(err_h - l2_err_h[i]) < 10.0^-12
  end
end

end # module

Executing @report_opt Main.Williamsom2ThetaMethodFullNewtonTests.test() then generated the attached file.

JET_report_opt_Williamson2ThetaMethodFullNewtonTests.txt

Highly related to gridap/Gridap.jl#646 and gridap/Gridap.jl#641

The current approach that we are using for generating the cubed sphere mesh has a bottleneck that will prevent us from scaling beyond hundreds of cells per panel (indeed, my laptop runs out of memory for approx 200x200 cells/panel).

The issue is in this line:

https://github.com/gridap/Gridap.jl/blob/a52bc10503f77f700b7cee096c7311aebc131cee/src/Geometry/CartesianDiscreteModels.jl#L53

Recall that we generate a Cartesian volume mesh of the cube, which we then restrict to the cube surface. In this particular line, the Cartesian Grid is (temporarily) transformed into a UnstructuredGrid (coordinates vertices and connectivity explicitly stored, instead of managed implicitly), and this causes memory to increase cubically with the number of cells per direction (i.e., if we multiply the number of cells per direction by 2, the memory grows by a factor of 8). This is clearly not acceptable, it should increase quadratically, as we have a 2D mesh.

Possible approaches that comes into my mind (perhaps there are better solutions):

• To bypass this issue is to write an ad-hoc extension of DiscreteModel{2,3} for the cubed sphere mesh made of 6 nxn CartesianDiscreteModel{2,2} panels joined together (as suggested by @davelee2804 in one of his comments) + the alpha/beta equiangular coordinates to R^3 coordinates that @cbladwell already tested previously.
• Use p4est to handle the cubed sphere mesh. This would be required for large scale simulations.

Hi @davelee2804 @cbladwell,

find below a code snippet to assemble the (possibly block diagonal) L2 mass matrix and its inverse "in one shot". I think this is an efficient way to go. Note that:

1. The code avoids extracting the diagonal blocks from an already assembled global sparse matrix. This is very inefficient. (scalar sparse matrix storage formats are not well prepared for such kind of operations).
2. The code collects the contribution of each cell to the global matrix into a regular vector (i.e., no longer lazy) using collect, i.e., explicitly storing all cell matrices at once (i.e., the matrix diagonal blocks). Usually we cannot afford this. However, in this particular case, the amount of memory required is equivalent to the one storing all global matrix entries. As a result, we avoid computing the cell matrices twice (once for L2MM and another one for invL2MM).
3. We use similar in order to get a sparse matrix with the same sparsity pattern as L2MM but with undefined entries yet. This is far more efficient than assembling a matrix from scratch. (We have already commented this in other contexts).

Let me know if everything makes sense ....

using Gridap
using GridapGeosciences

order=1
degree=4

model = CubedSphereDiscreteModel(1; radius=rₑ)

Ω = Triangulation(model)
dΩ = Measure(Ω, degree)

reffe_lgn = ReferenceFE(lagrangian, Float64, order)
Q = FESpace(model, reffe_lgn; conformity=:L2)
P = TrialFESpace(Q)

function assemble_L2MM_invL2MM(U,V,dΩ)
  a(a,b) = ∫(a⋅b)dΩ
  v = get_fe_basis(V)
  u = get_trial_fe_basis(U)
  assemblytuple    = Gridap.FESpaces.collect_cell_matrix(U,V,a(u,v))
  cell_matrix_MM   = collect(assemblytuple[1][1]) # This result is no longer a LazyArray
  newassemblytuple = ([cell_matrix_MM],assemblytuple[2],assemblytuple[3])
  a=SparseMatrixAssembler(U,V)
  L2MM=assemble_matrix(a,newassemblytuple)
  for i in eachindex(cell_matrix_MM)
     cell_matrix_MM[i]=inv(cell_matrix_MM[i])
  end
  invL2MM=similar(L2MM)
  Gridap.FESpaces.assemble_matrix!(invL2MM, a, newassemblytuple)
  L2MM, invL2MM
end

A,B=assemble_L2MM_invL2MM(P,Q,dΩ)
println(A*B)