clever steady state

We can formulate the equations of motion as ∂tu+F(u)=f(u) \partial_t \mathbf{u} + \mathbf{F}(u) = \mathbf{f}(u)

We wish to find a steady state, namely ∂tu=0\partial_t u = 0 which gives us the relation F(u)−f(u)=0 \mathbf{F}(u) - \mathbf{f}(u) = 0

Let G(u)=F(u)−f(u) \mathbf{G}(u) = \mathbf{F}(u) - \mathbf{f}(u)

Then we can root find using newton's method, namely un+1=un−JG(u)−1G(u) \mathbf{u}_{n+1} = \mathbf{u}_n - J_\mathbf{G}(\mathbf{u})^{-1} \mathbf{G}(u)

Here's the difficulty: how do we compute JGJ_{\mathbf{G}} ? The beauty of using a spectral method (c.f. the Guerra and Ullrich paper) is that the method is global and this Jacobian is assumed to be dense. This induces a substantial cost to increasing the order of your method; however, the computing derivatives with respect to spectral coefficients is downright trivial. I much prefer the parallel scaling of the spectral element dycore, however I'm unsure how computing such derivatives interacts with the DSS (so-called direct stiffness summation) step that's necessary to project solutions back into continuous Galerkin space. However, the HOMME codebase appears to have fully implicit solver infrastructure, which means that there may be a template for how to do this in a principled way.

Solving the linear equation is akin to any elliptic (I think) lumped FEM method which means it would be sensible to do something like preconditioned GMRES.

Using clues from HOMME:

In the file ${HOMME_ROOT}/share/prim_implicit_mod.F90 we find a clue in the subroutine residual, which demonstrates the infrastructure to solve F(u)=(un+1−un)Δt−DSS [RHS(un+1)] \mathbf{F}(\mathbf{u}) = \frac{(\mathbf{u}_{n+1} - \mathbf{u}_{n})}{\Delta t} - \textrm{DSS }\left[ \textrm{RHS}(\mathbf{u}_{n+1}) \right]

Under our assumptions we would revise this equation so un+1=un\mathbf{u}_{n+1} = \mathbf{u}_{n}, i.e.

F(u)−DSS [RHS(un+1)]=0 \mathbf{F}(\mathbf{u}) - \textrm{DSS }\left[ \textrm{RHS}(\mathbf{u}_{n+1}) \right] = 0

which seems like it may be a straightforward modification to make to the code?

I also note that in the C++ utilities, the Nox trilinos library is included, which is a general purpose nonlinear system solver that works with Kokkos.

NOTE! ${HOMME_ROOT}/src/implicit_mod.F90 contains a bunch of interface routines to Nox! This should work swimmingly.

The key is understanding noxsolve, which appears to only be included in the preqx codebase.

The picard linearization subroutines are where it's at! possibly useful stuff about picard linearization

Note: remove all dependence on time AND remove any mention of hyperviscosity.

Modifying subroutines

    void (*precFunctionblock11)(double*, int, double*, void*)      = sw_picard_block_11;
    void (*precFunctionblock12)(double*, int, double*, int, void*) = sw_picard_DFinvBt;
    void (*precFunctionblock21)(double*, int, double*, int, void*) = sw_picard_block_21;
    void (*precFunctionblock22)(double*, int, double*, void *)     = sw_picard_schur;

    void (*get_globalIDs)(int, int *, void *) = homme_globalIDs;
    void (*get_HelmElementMat)(int, int, double *,int *, void *)=helm_mat;
    void (*get_HelmMap)(int, int, int *, void *)=helm_map;

Another quantity of interest: up is the u at the next timestep.

IF I want to try this with theta-l, I can use the ARKODE and sub out KINSOL to do what I want to do

Case configuration

Next up! Setup a sanity check case. Set u = 0 to start.

I.e. specify u=0,v=0,T=T0,ps=p0,zs=0m u=0, v=0, T=T_0, p_s = p_0, z_s = 0\mathrm{m}

Pick a DCMIP test to bastardize