slim3D: add export netCDF + full + restart

dgpy/scripts/slim3d.py

@@ -158,18 +158,20 @@ class Slim3d_equations :
 
 class Loop:
     """Temporal solver"""
-    def __init__(self, equations, maximum_2d_time_step=3600, ratio_3dvs2d_time_step=1, export_time_step=3600, initial_time='1970-01-01 00:00:00', final_time=0, output_directory="./output"): 
+    def __init__(self, equations, maximum_2d_time_step=3600, ratio_3dvs2d_time_step=1, export_time_step=3600, ratio_full_export=-1, initial_time='1970-01-01 00:00:00', final_time=0, output_directory="./output", time_zone=0): 
         self._slimSolver = equations._slimSolver
         self._equations = equations
         self._max_2d_dt = maximum_2d_time_step
         self._dt2d3d_ratio = ratio_3dvs2d_time_step
         self._export_dt = export_time_step
+        self._ratio_full_export = ratio_full_export
         fmt = '%Y-%m-%d %H:%M:%S'   
         date = slim_private.datetime.datetime.strptime(initial_time, fmt)
-        self._initial_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
+        self._initial_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second]) + time_zone * 3600.
+        self._initial_date = date
         if slim_private._is_string(final_time):
             date = slim_private.datetime.datetime.strptime(final_time, fmt)
-            self._final_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
+            self._final_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second]) + time_zone * 3600.
         else:
             self._final_time = self._initial_time + final_time
         self._odir = output_directory
@@ -181,8 +183,16 @@ class Loop:
         self._export_eta      = False
         self._export_kappav   = False
         self._export_nuv      = False
+        self._export_uv_nc    = False
+        self._export_S_nc     = False
+        self._export_eta_nc   = False
+        self._export_z_nc     = False
         self._timeIntegrator = None
         self._export = []
+        self._ncExport = []
+        self._restart_dir = None
+        self._restart_ind = -1
+        self._time_zone = time_zone
     
     def export_uv(self): self._export_uv = True
     def export_w(self): self._export_w = True
@@ -192,19 +202,38 @@ class Loop:
     def export_elevation(self): self._export_eta = True
     def export_vertical_diffusivity(self): self._export_kappav = True
     def export_vertical_viscosity(self): self._export_nuv = True
+    def export_uv_nc(self): self._export_uv_nc = True
+    def export_salinity_nc(self): self._export_S_nc = True
+    def export_elevation_nc(self): self._export_eta_nc = True
+    def export_z_nc(self): self._export_z_nc = True
+        
+    def restart(self, directory, index):
+        self._restart_dir = directory
+        self._restart_ind = index
 
     def setup(self):
         slim_private.slim3d_setup(self)
         slimSolver = self._slimSolver
         if dgpy.Msg.GetCommRank() == 0:
-            if not slim_private.os.path.exists(self._odir):
-                slim_private.os.mkdir(self._odir)
+            if not slim_private.path.exists(self._odir):
+                slim_private.makedirs(self._odir)
         dgpy.Msg.Barrier()
         d = slimSolver.getDofs()
+        self._ncWriter2d = dgpy.slimNetCDFIO(slimSolver.groups2d)
+        self._ncWriter2d.setInitTime(self._initial_date.year,self._initial_date.month,self._initial_date.day,self._initial_date.hour,self._initial_date.minute,self._initial_date.second,self._time_zone,00)
+        self._ncWriter3d = dgpy.slim3dNetCDFIO(slimSolver.columnInfo,self._ncWriter2d)
+        if self._export_eta_nc: self._ncExport.append([d.etaDof2d, self._odir + '/eta','eta'])
+        if self._export_uv_nc: self._ncExport.append([d.uvDof, self._odir + '/uv','uv'])
+        if self._export_z_nc: self._ncExport.append([d.zDof, self._odir + '/z','z'])
         if self._export_uv:     self._export.append(dgpy.dgIdxExporter(d.uvDof, self._odir + "/uv", True))
         if self._export_w:      self._export.append(dgpy.dgIdxExporter(d.wDof3d, self._odir + "/w"))
         if self._export_uvAv2d: self._export.append(dgpy.dgIdxExporter(d.uvAvDof2d, self._odir + "/uvAv2d", True))
         if self._export_eta:    self._export.append(dgpy.dgIdxExporter(d.etaDof2d, self._odir + "/eta"))
+        if self._export_S_nc: 
+            if slimSolver.getSolveS():
+                self._ncExport.append([d.SDof, self._odir + '/S','S']) 
+            else: 
+                dgpy.Msg.Warning("Salinity is not solved. It won't be exported")       
         if self._export_S: 
             if slimSolver.getSolveS():
                 self._export.append(dgpy.dgIdxExporter(d.SDof, self._odir + "/salinity"))
@@ -225,17 +254,22 @@ class Loop:
                 self._export.append(dgpy.dgIdxExporter(d.nuvDof, self._odir + "/nuv"))
             else: 
                 dgpy.Msg.Warning("Turbulence is not solved with gotm. nu_v won't be exported")
+        return self._timeIntegrator
 
     def loop(self):
-        slimSolver = self._slimSolver
-        d = slimSolver.getDofs()
         if not self._timeIntegrator:
             self._timeIntegrator = self.setup()
         timeIntegrator = self._timeIntegrator
-
+        slimSolver = self._slimSolver
+        d = slimSolver.getDofs()
+        if self._restart_ind > -1:
+            timeIntegrator.reStart(self._restart_dir, self._restart_ind)
         if dgpy.Msg.GetCommRank()==0 : dgpy.Msg.Info("export %i" % timeIntegrator.getExportCount())
         for e in self._export :
             e.exportIdx(timeIntegrator.getExportCount(), timeIntegrator.getTime() - self._initial_time)
+        for dof,fname,field in self._ncExport :
+            filename = fname+'_{0:05d}'.format(timeIntegrator.getExportCount())
+            self._ncWriter3d.exportUGRID(dof, filename, timeIntegrator.getTime() - self._initial_time, timeIntegrator.getIter(),field)
         tic = slim_private.time.time()
         while timeIntegrator.getTime() < timeIntegrator.getEndTime() :
             timeIntegrator.advanceOneTimeStep()
@@ -244,6 +278,11 @@ class Loop:
                 if dgpy.Msg.GetCommRank()==0 : dgpy.Msg.Info("export %i" % timeIntegrator.getExportCount())
                 for e in self._export :
                     e.exportIdx(timeIntegrator.getExportCount(), timeIntegrator.getTime() - self._initial_time)
+                for dof,fname,field in self._ncExport :
+                    filename = fname+'_{0:05d}'.format(timeIntegrator.getExportCount())
+                    self._ncWriter3d.exportUGRID(dof, filename, timeIntegrator.getTime() - self._initial_time, timeIntegrator.getIter(),field)
+            if timeIntegrator.checkFullExport() :
+                timeIntegrator.exportAll()
             timeIntegrator.checkTracerConsistency()
             uvNorm = d.uvDof.norm()
             if dgpy.Msg.GetCommRank() == 0 : dgpy.Msg.Info("%5i [ %s ] normuv %6.12e clock %.1fs" % (timeIntegrator.getIter(),slim_private.time.strftime("%d %b %Y %H:%M:%S", slim_private.time.gmtime(timeIntegrator.getTime())), uvNorm, slim_private.time.time() - tic))

dgpy/scripts/slim3d_setup.py

@@ -328,8 +328,10 @@ def slim3d_setup(loop):
     dt = CFL*dtRK
     maxUV = 1.5 # max advective velocity
     c2d = slim_private.math.sqrt(9.81*20.)
-    dt = c2d/(c2d+maxUV) * dt
-    dt = min(dt, loop._max_2d_dt)
+    dt1 = c2d/(c2d+maxUV) * dt
+    dt2 = min(dt1, loop._max_2d_dt)
+    n_iter_per_export = numpy.ceil(loop._export_dt/(dt2*loop._dt2d3d_ratio))
+    dt = loop._export_dt/(n_iter_per_export*loop._dt2d3d_ratio)
     if dgpy.Msg.GetCommRank() == 0 :
         dgpy.Msg.Info('Runge-Kutta dt :' + str(dtRK) + ' CFL:' + str(CFL) + ' dt:' + str(dt))
     
@@ -341,7 +343,7 @@ def slim3d_setup(loop):
     timeIntegrator.setExportInterval(loop._export_dt)
     #timeIntegrator.setExportAtEveryTimeStep(False)
     #timeIntegrator.setFullExportAtEveryTimeStep(False)
-    #timeIntegrator.setFullExportIntervalRatio(100)
+    timeIntegrator.setFullExportIntervalRatio(loop._ratio_full_export )
     timeIntegrator.setTimeStep(dt)
     timeIntegrator.set3dTimeStepRatio(loop._dt2d3d_ratio)
     timeIntegrator.setInitTime(loop._initial_time)