Merge pull request #1 from colinjcotter/SSPRK3

Version used for QJRMS 2023 paper.

README.txt

@@ -0,0 +1,6 @@
+The version used to produce results for the qjrms paper.
+
+averaged_sw_explicit.py: the averaged model
+standard_sw.py: the standard model
+calculate_norm.py: compare the .h5 files produced by the above models and calculate norms
+qjrms/ directory: contains data files and scripts which were used to produce figures 4-6 of the paper

averaged_sw.py

@@ -58,15 +58,17 @@
 R0 = 6371220.
 H = Constant(5960.)
 
+mesh_degree = 3
 mesh = IcosahedralSphereMesh(radius=R0,
-                             refinement_level=ref_level, degree=3,
+                             refinement_level=ref_level, degree=mesh_degree,
                              comm = ensemble.comm)
 cx = SpatialCoordinate(mesh)
 mesh.init_cell_orientations(cx)
 
 cx, cy, cz = SpatialCoordinate(mesh)
 
-outward_normals = CellNormal(mesh)
+outward_normals = interpolate(CellNormal(mesh),VectorFunctionSpace(mesh,"DG",mesh_degree))
+
 perp = lambda u: cross(outward_normals, u)
 
 V1 = FunctionSpace(mesh, "BDM", 2)
@@ -199,7 +201,7 @@
 lambda_c = -pi/2.0
 phi_x = asin(x[2]/R0)
 phi_c = pi/6.0
-minarg = Min(pow(rl, 2), pow(phi_x - phi_c, 2) + pow(lambda_x - lambda_c, 2))
+minarg = min_value(pow(rl, 2), pow(phi_x - phi_c, 2) + pow(lambda_x - lambda_c, 2))
 bexpr = 2000.0*(1 - sqrt(minarg)/rl)
 b.interpolate(bexpr)
 
@@ -211,7 +213,7 @@
 DU = Function(W)
 V = Function(W)
 
-U_u, U_eta = U.split()
+U_u, U_eta = U.subfunctions
 U_u.assign(un)
 U_eta.assign(etan)