Referee statement

I suggest to publish this article in GMD journal after Minor revisions.

General remarks

The manuscript is well written and also well constructed, it is then easy to follow the various cases

tested by the authors, from simple 2D objects to more complex objects. The main objective of this

manuscript is to demonstrate Firedrake’s applicability for geodynamical simulation. The first part of

the manuscript focus of reproducing already approved analytical test cases and well known benchmarks.

The authors must also be congratulate for the efforts they have put in describing the previous work

done on the topic, which later help understanding the conducted experiments and results.

The choice made by the authors to present their Python code by section, displaying each time only

differences compared to the initial test case must be acknowledge, as it is very helpful in understanding

the paper, and a great asset in reproducibility.

Major comments

Although the paper is well constructed and authors guideline conducting their test cases is relatively

easy to understand, I believe the paper readability could be enhanced with the addition of a global

table introducing the various cases(from the 2-D to a realistic application in 3-D spherical geometry),

their differences and their interests.

Some improvement could be brought to figures, in many cases the Y axes could homogeneous to ease

reader analysis of the results (cf. Figure 1, 7, 8...).

Minor comments

None.

I have read with great interest the paper by Davies et al on Firedrake in geodynamics. I found the paper well written and well structured. I also really appreciate that the authors explain how the codes work almost line by line while presenting results of standard benchmarks in the community. 

I have a few comments:

- section 3 is called 'Governing equations'. As such I would expect that the generic form of the Stokes equations and the energy equation are then introduced. However equations 1&2 are dimensionless and those of an incompressible fluid under Boussinesq approximation. These are of course the relevant equations and the necessary ones for the Blankenbach et al benchmark, but maybe Section 3 should start with a foreword that the paper only deals with buoyancy driven flow of a single fluid so that the relevant set of equations are eqs 1&2.  

- Given how firedrake works, i.e. the user must input the weak form of the equations directly in the code, I appreciate that the authors have decided to quickly re-derive the weak form of the equations because these find their way in Listing 1. However I am rather puzzled by their approach for the weak form of the continuity equation (eq. 9). Technically this equation is correct, especially when Eq. 13 is taken into account, but no justification is given as to why div(u)=0 is integrated by parts. This is highly uncommon, so much so that I was not able to find a reference/textbook in which this is carried out. At the very least I believe the authors should justify why this uncommon approach is taken. Conversely, would the 'regular' way of not integrating the continuity equation by parts also work in Firedrake? Also, if this form of the weak form (eq 9) is used this begs the question as to what the spaces U,W are but these are not specified beyond "we select appropriate function spaces" at line 203.

- At line 211, I am not sure to understand what the authors mean by stating that the weak forms are "a more general mathematically *rigorous* representation of the continuous PDEs". In what sense?

- in section 4.2 the authors state "To simplify we will solve for velocity and pressure, u and p, in a separate step before solving for the new temperature T n+1." It migh be beneficial in the case to recall first that Eqs. 8-10 are coupled and how, and that one can solve Stokes followed by the energy equations but that this process should be iterated out. 

- line 490-491: equation 36 is not formulated as the second invariant of strain rate since the second moment invariant is \frac12 epsilon:epsilon . The denominator is then the square root of twice the moment invariant in this case. May be calling it effective strain rate could be better?

- line 504: MUMPS is always capitalised in the MUMPS literature

- line 524, line 625. For information, this setup is already implemented in ASPECT (see manual). You could use ASPECT results as a point of comparison too.

- section 5.3.1: although it does not matter for the steady state, there is no mention of the erf function in Case 1 of Busse et al 1994. Also, maybe already mention that advection would need a stabilisation term for high Peclet numbers (which is discussed later) but in this specific case it is not needed? finally, this problem only contains a single pressure nullspace since Dirichlet bc are prescribed all around, so I am not following point 5 at line 580, but it probably comes from my lack of understanding what a near-nullspace is. May be worth quickly explain what these are (and how they differ from nullspace) in section 4.3 and why they matter so much (effectively improving the last paragraph of the section by a sentence or two). Could it also be an idea to run one model without this near-nullspaces information passed to the GAMG preconditioner and show how much it matters ? (in an appendix?) 

- Section 5.3.3: This benchmark was run with ASPECT in Liu & King, doi:10.1093/gji/ggz036

- general comment about Nu and vrms calculations throughout the paper: may be show/discuss how these are computed and how it translates into firedrake code ?

 

- fig 10: I would have expected the scale on the x-axis of these plots to be logarithmic

-line 723: the 'trilinear' adjective alone can't be used to describe 'Q2-Q1' 

- Eq.40: I find that this is a bit of a missed opportunity. Given the type of modelling in Section7, why not use the radial viscosity profile of Bull et al (EPSL 2014,  http://dx.doi.org/10.1016/j.epsl.2014.06.008) who also drive the system with plate velocities and an identical temperature dependence and briefly compare/discuss your results and theirs ? (although I understand that the focus of the paper is firedrake, not plate tectonics)

- Section 7: I think it would be nice to report on how long that simulation took on how many cores. 

After reading the paper I am impressed and convinced that Firedrake has a role to play in computational geodynamics in the coming years. However, the numerical experiments here, despite being useful as benchmarks and learning cases, are limited: they do not showcase features which have now become very common in state-of-the-art such as free surface, elasto-visco-plastic rhelogies, particle-in-cell or level set or compositional fields, multiphysics with operator splitting, etc ... I think the authors should address this point in the discussion by explaining whether they think that Firedrake as it is now is capable to generate a code with the above features, what would/could be the problems along the way, etc ... For instance the use of compositional fields requires to solve an advection equation per field, which obviously Firedrake could do. But what about particle-in-cell ? Which limitations (if any) will Firedrake impose on what can and cannot be implemented?

Finally, on the basis of the benchmarks presented I would argue that the authors have demonstrated that Firedrake is suitable for single fluid convection modelling, not for lithospheric/crustal modeling. May be the title should reflect this and the word 'geodynamics' be replaced by Mantle convection?

Could the authors also add an appendix with vrms, Nu values in tables ?

Cedric Thieulot.

I've read the paper about using Firedrake for mantle convection with

great interest. It is clear that Firedrake is an excellent system for

writing such simulators in quite a small amount of code, something I

think is just really impressive!

That said, there are two big issues I would like to raise (and on

which I will then let the editor make a decision):

* What is it actually that is new/what is the *message* of the paper?

* The paper is quite long.

Let me take these on in turn:

1/ It hasn't quite become clear to me what the *message* of the paper

is or should be. Most of the paper is devoted to things that are quite

standard for those who have used the finite element method, and to

showing numerical results. But neither of that is new: We know how the

finite element method works. And how we turn the weak form into a

linear system, and one would expect that the authors have checked that

their code produces correct answers -- the fact that it does could

have been stated in a couple of sentences of the form "We have

verified that our code produces correct answers for benchmarks X, Y,

and Z, to within the errors that have previously been reported when

using other codes (ref 1, ref 2, ...)". We don't need to see many

pages of tables and figures to believe this; this is particularly true

because these figures do not actually compare against other codes, and

so all we can infer from them is that the code is correct, but not

whether it is better or worse in any regard than existing codes.

So, then, one is left with the question what actually *is* the message

of the paper, seeing that the authors do not actually use the space

well to make their point. I *suspect* that the authors think that the

message should be "Using Firedrake, we can write these codes in 200

lines that using other software systems takes 2000 or 20,000". If this

were the message, I think it would make a fine paper, though then I

would try to remove much of the extraneous stuff mentioned above (see

also my point 2 below). 

I would, however, like to point out that that by itself is not as

impressive as it may sound like. Being able to do what others have

been able to do for a long time is not actually particularly

interesting, even if the new system can do it quicker and with less

code. The question that *should* have been answered is what one can do

with this system that others cannot do -- what does Firedrake make

possible that would otherwise not have been possible (including

because it would take too much work in other systems)? Unfortunately,

the paper does not answer this question, and I'm not sure whether

there really *is* a good answer. The fundamental issue is that systems

like Firedrake are really good at *making simple things easy*; the

paper illustrates this: all of the codes shown are very nice and

concise, but they fundamentally all solve very simple problems that

the community has been able to solve for a long time. But the other

codes that are out there can *also* solve much more complicated

things, with moving boundaries, much more complex material laws,

adaptive meshes, adaptive time step choices, and particles. It is not

a stretch to speculate that doing all of these things would not

actually be much easier with Firedrake either -- the majority of the

code in ASPECT, for example, is after all not in the description of

the finite element shape functions, or the description of linear

solvers, but in material models, postprocessors, time step control,

dealing with how compositional fields affect this that and the other,

etc. It wouldn't surprise me if that could be done with 5x less code

in Firedrake than in C++, but that would still take several 10,000

lines of code. (I will add as a sidenote that the model of writing

many small codes for each benchmark is also not sustainable in that it

encourages the fragmentation mentioned early on in the paper where

everyone has their own variation of the code for a specific problem;

codes like ASPECT spend many thousands of lines of code on run-time

parameter handling precisely so that we can have a single code base

and only need to exchange input files.)

In other words, I would appreciate if the authors could sharpen their

message: Being able to do what others have already done, and

illustrating this with a long list of examples is just not very

interesting in itself.

2/ Space: The paper is quite long, at 50 pages. That's because in many

regards it reads like a *manual*: It is addressed to people who don't

already have the finite element background to understand the basics

that underly the codes, and to people who need to be convinced with a

long list of detailed examples that the codes produce the correct

answer. But the audience for a manual is different from the audience

of a research paper. The authors might want to think about who they

are writing this paper for, and then think about every section and

every figure, and re-evaluate whether that section or figure is

necessary for that audience. For example, I thought that many of the

figures carry no more meaning than one could equally express with the

sentence "We have verified that our simulations converge to the same

values as reported in X and Y (1997)". It would not greatly surprise

me if the paper could be shrunk to 30 or fewer pages without actually

losing any of the message.

Other comments -- only a few, the paper is actually really well written:

l. 206: This is a rather unusual formulation, with not integrating by

parts the 'grad p' term, but integrating by parts the 'div u'

term. The only comparable formulation I know of is what is done when

comparing the primal and dual mixed formulations of the mixed Laplace

equation -- there it makes sense based on whether one wants the

velocity or the pressure to be in a space that requires

derivatives. But for the Stokes equations, the velocity *always* needs

derivatives because the viscous friction term already has derivatives

on it. The usual approach is to get the derivative off the pressure

and let it be in L2, but the authors here gratuitously require the

pressure to be in H1. This seems unnecessary and at least requires

explanation.

l. 277: Choosing C_ip=100 seems like a sledgehammer. It will for sure

affect the condition number of the matrix. Is there not a smarter

approach to choosing it.

l. 284: Choosing Crank-Nicolson isn't stupid, but it is also not

particularly good choice. Why not use something a bit more accurate?

Or, if you want to use something simple for expository reasons, at

least say so. Separately, all codes use a fixed time step size; this

too is not sufficient in practice.

l. 323: For more complex rheologies -- specifically if the viscosity

also depends on the pressure -- the Newton matrix will have additions

also in other blocks. But there are other complications one has to

address for nonlinear rheologies as well; a straight up Newton scheme

does not always work for strain-weakening rheologies as are common in

geodynamics (see Fraters et al., GJI, 2019). In fact, this falls in

the category of things that just require a decent amount of code: One

has to choose a Newton matrix that is *not* just the derivative of the

residual, and that has to be programmed because the

auto-differentiation of the residual isn't going to produce the matrix.

p. 13: For this and all other listings: Align comments vertically to

make the code easier to read. If you can align a few '=' signs, then

that's worth doing as well.

l. 497: "of of" -> "of"

l. 683, and the following paragraph about solver scaling: This does

not actually sound true. It is simply a deficiency of the

preconditioner the authors choose. There are plenty of papers that

show that with good preconditioners, one can achieve almost perfect

weak scaling, and in particular achieve a more or less constant number

of linear iterations regardless of problem size. The statement as

given is a cop-out :-)

Fig. 10: Please use a logarithmic x-axis to make these graphs easier

to read.

Also Fig. 10: It would have been nice if there was a comparison to how

long these sorts of computations take with other software systems. For

example, how long would ten time steps with ASPECT take? This goes in

the same direction as the "mission" of the paper I mentioned

above. Just knowing *that* your code can solve a problem others have

been able to solve for a long time is not so interesting. If you can

do it *as fast as others* with 1/100th the lines of code they need,

then that *is* interesting.

Section 7 (l. 720): It doesn't sound quite right to call this section

"realistic applications". The application uses an incompressible

mantle, a linear temperature and depth dependence for the viscosity,

and has no special provisions to deal with the crust any different

than the rest of the mantle. This is neither realistic, nor

new. "Realistic" models have been far more complex than this for a

very long time already.

l. 733: Are T=0 and T=1 mixed up here?

This is a well-written manuscript. I would recommend it to be published with some minor revisions.

Major issues:

1) I agree with RC4 by Wolfang Bangerth that the main message of the manuscript should be clarified. If the target audience is geodynamicists who are testing different approximations of governing equations or the performance of different solver algorithms, this manuscript is an excellent introduction and demonstration and can be published with some minor revisions in the introduction and conclusion. If the target audience is geodynamicsists who are working on real-world geophysical problems, this manuscript falls short in several aspects. It is probably beyond the scope of this manuscript to address these aspects, but I would list them below nonetheless:

a) Since solving the governing equation is only a small part of a mantle convection code, some part of the generated C kernel must be extended. I would like to know how easy is it to modify and extend the generated C kernel? Can it be done in python or must be done in C? For example, to add dynamic time stepping, one has to calculate min(v/h) over all elements, where h is the element size. Can the contents of v and h be accessed in python, or only in C? 

b) Most legacy mantle convection codes are using domain decomposition to parallelize. But Firedrake seems to parallelize by distributing the matrixes and vectors. How will this affect the extension of the code, for example, adding markers?

Some minor suggestions:

2) In the code listings, it will be better if some function arguments are passed as keywords, instead of just as numbers. For example, line 8 in listing 1 becomes: W = FunctionSpace(mesh, family="CG", degree=1)

3) Line 40-41 linsting 1, the boundary tags 1, 2, 3 ,4 are better to named as left_id, right_id, bottom_id, top_id, respectively. Related, Listing 2, bottom_id and top_id are used but not defined.

4) The solver parameters for energy solver in Listing 1 and as described in line 405 appear to be inconsistent with those described in line 300. Also, the solver parameters for Stokes solver in line 336 are inconsistent with those in later code listings.

5) In the code, the buoyancy term is using Ttheta, rather than Told or Tnew. Why?

I agree with the comments RC4 (by Wolfgang Bangerth) and RC5 that the overall purpose of the paper needs to be clarified.

In general the authors should try to answer the following question: What can we do with this code now that was impossible with other approaches before?

I'll try to provide some suggestions for the authors in order to answer this questions and improve the manuscript to a point where I can recommend it for publication.

Major remarks:

The authors actually give some suggestions themselves in Section 8 of the paper which have not been followed up but would have made the manuscript

a clear candidate for publication. Furhtermore, I would suggest to restructure the paper such that we have a section presenting the benchmark implementations using UFL and a part discussing improvements over other existing packages and approaches.

Benchmark cases implemented using UFL: It would be a valuable contribution to provide a ready-to-use collection of benchmark cases (as presented in the paper) described in UFL to be used with multiple packages that can make use of UFL.  An example here is dolfin_dg (https://bitbucket.org/nate-sime/dolfin_dg/),  which provides UFL forms for discontinuous Galerkin discretizations for compressible flow and can be used by Fenics, Firedrake, or even DUNE-FEM. Right now it is not entirely clear whether the presented examples could be used with Fenics or TerraFERMA. It would also be in line with the much advertized separation of concerns.

Improvements over other existing approaches: In the second part of the paper I suggest to better highlight the strength of using the benchmark cases with Firedrake over other available packages. For example:

• Could we easily define our own preconditioner, for example, making use of UFL again? And if so, an example should be added.
• The presented finite element discretization (Q2-Q1 on Cartesian grids) is standard. How easy and feasible is it to use other (maybe more appropriate) discretiztions or other grid element types?
• Could we easily do a space-time discretization since UFL should make it actually very easy to write this down?
• As stated by the authors, "the automated approach underpinning Firedrake has the potential to revolutionize the use of adjoints and other inverse schemes in geodynamics". This would have been a clear major improvement over the existing state of the art. Could this at least be demonstrated for a simplified example? It would clearly help the community to have a starting point in this direction.
• Scalability: Would it be possible to also see at least one strong scaling result for one of the presented cases?

Limitations: What are the limitations of the presented approach? Existing approaches already cover a wide range of features.For example, could this code be used to replicate the results presented in "Large-scale adaptive mantle convection simulation" by Burstedde, Stadler et al. (https://doi.org/10.1093/GJI%2FGGS070) at the same grid resolution?

Besides that, the length of the paper in it's current form is ok, since for this journal the resulting published paper will be much shorter due to the two-column format used. I'm just wondering how the code representation would look like in that format. And then, addressing my comments from above will probably mean to shorten the current presentation of the test cases. In addition, I have a few minor remarks.

Minor remarks:

1. Throughout the paper: Crank-Nicholson should be Crank-Nicolson.
2.  Line 418: The citation of Homolya et al. is in my opinion wrong here, because the cited paper in the end only shows compilation times of the TSFC but does not make any statement about the efficiency of the produced code.

This is a comprehensive and detailed paper on implementing mantle convection simulations in Firedrake. I am in favor of publishing it, due to the attention to reproducibility, the careful matching of benchmarks, and its general model and tutorial value. I have only minor suggestions:

The Newton terms are omitted, the Nitsche terms are omitted form the discussion of system solution. Yet these would be rather interesting since these take some effort to derive and many will not have taken the trouble to do so before, especially in light of concrete benchmarks. Conversely, what is presented in terms of equations in the paper is rather well known. Thus my question, would it be practical to add these equtations/derivations? Or is it rather that due to the automatic differentiation available in Python, there is no need to present them in detail?

Figure 6b is the first that disagrees with established values.  Is it understood why?  Maybe the results obtained here are actually better?

Chapter 7, realistic convenction simulation: what is the spatial and temporal resolution?

Line 220, it is also required that finite element of different spaces pairings are stable, which depends crucially on the choice of these spaces. The paper omits this point entirely, maybe some comments can be added. Related, discontinuous pressure discretizations may be worth attention.

Line 255: The (imposed) Tinhom values produce a term linear in q that moves to the right hand side of (10).  How is this handled?

Line 290, 298: It is not clearly described how the temperature equation could possibly become nonlinear in T.  Maybe the authors are referring to the fully coupled systems, where the parameter u depends itself on T?

Line 550: Using petsc fieldsplit is more advanced and even more interesting than what is lined out earlier in the paper.  It would be valuable to discuss the corresponding numerical approach and technique in more detail, which may help readers to move in this direction themselves.

Line 836, may mention precondition work by Bangerth used in the Aspect code.