ekman

benchmarks/ekman_prod/test.py

+import dgpy
 import slim3d
 import slimPre
 import numpy
 import sys
 
 Ek = 0.05
-nu = 0.01
-f = 1e-4
-H = (nu/(f*Ek))**0.5
 D = 1e-5
 C = 1e-5
-nb_layers = 10
+
 if len(sys.argv) == 4 :
     Ek = float(sys.argv[1])
     C = float(sys.argv[2])
     D = float(sys.argv[3])
 
-print("H %.16g" % H)
-print("Ek %.16g" % Ek)
-print("nu %.16g" % nu)
-print("f %.16g" % f)
-print("D %.16g" % D)
-print("C %.16g" % C)
-print("nb_layers %d" % nb_layers)
+nu = 0.01
+f = 1e-4
+H = (nu/(f*Ek))**0.5
+nb_layers = 10
+
 
 
 ### Prepro ###
@@ -42,7 +38,6 @@ taux = (D/2*xyz[:,0]-C/2*xyz[:,1])
 tauy = (C/2*xyz[:,0]+D/2*xyz[:,1])
 slimPre.write_file('windstress.nc', region=region_global, time=None, data=[('taux',taux),('tauy',tauy)])
 
-
 ### Run ###
 domain = slim3d.Domain('box.msh', reference_density=1027)
 equations = slim3d.Slim3d_equations(domain, temperature=False, salinity=False)
@@ -64,20 +59,28 @@ equations.set_horizontal_viscosity(mode="smagorinsky")
 
 time_loop = slim3d.Loop(equations, 
         time_step=400,
-        export_time_step=0.5*3600,
+        export_time_step=3600,
         final_time=200*24*3600, 
         output_directory='output')
 
 time_loop.export_elevation()
-time_loop.export_uv(True)
-time_loop.export_w()
+time_loop.export_uv(False)
+#time_loop.export_w()
 time_loop.export_uv2d(True)
 
 
 ################
 # vorticity 2d #
 ################
-import dgpy
+vortf=open("vorticity_log","w")
+vortf.write("H %.16g\n" % H)
+vortf.write("Ek %.16g\n" % Ek)
+vortf.write("nu %.16g\n" % nu)
+vortf.write("f %.16g\n" % f)
+vortf.write("D %.16g\n" % D)
+vortf.write("C %.16g\n" % C)
+vortf.write("nb_layers %d\n" % nb_layers)
+vortf.close()
 groups2d = dgpy.dgGroupCollection("box.msh", 1, None, 2, ["top_In","top_Bnd"])
 uvdof2d = dgpy.dgDofContainer(groups2d, 2)
 etadof = dgpy.dgDofContainer(groups2d,1)
@@ -87,7 +90,6 @@ dof = dgpy.dgDofManager.newCG(groups2d,2,sys)
 l2p = dgpy.L2ProjectionContinuous(dof)
 ucg = dgpy.dgDofContainer(groups2d,2)
 
-
 etagradf = dgpy.dgFEGradient(etadof.getFunction())
 etagradf.setBoundarySymmetry("Coast")
 etagrad = dgpy.dgDofContainer(groups2d,3)
@@ -137,6 +139,10 @@ def compute_vorticity(i) :
 
     print("2d  ",i,i2,j2)
     print("geo ",i,ig,jg)
+    vortf=open("vorticity_log","a")
+    vortf.write("2d  %i %g %g\n"%(i,i2,j2))
+    vortf.write("geo %i %g %g\n"%(i,ig,jg))
+    vortf.close()
     return i2
 ################
 
@@ -150,9 +156,13 @@ while time_loop.get_time() < time_loop.final_time:
         time_loop.print_iter_info()
         time_loop.export_fields()
         inew = compute_vorticity(iout)
-        if abs(iold-inew) < abs(inew)/100 :
+        if abs(iold-inew) < abs(inew)/1000 :
             break
         iold = inew
         iout += 1
+print("sucess")
+vortf=open("vorticity_log","a")
+vortf.write("success\n")
+vortf.close()
 time_loop.terminate()