slim3d: new time integrator

CMakeLists.txt

@@ -18,6 +18,7 @@ option(ENABLE_PETSC4PY "Enable PETSc4Py" ON)
 option(ENABLE_FUNCTIONC "ENABLE C CALLBACK" ON)
 
 project(dg CXX)
+set( CMAKE_EXPORT_COMPILE_COMMANDS 1 )
 set(CMAKE_CXX_STANDARD 11)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)

benchmarks/swe3d/tanganyika/2dxz/run.py

+import slim3d
+from slimPre import make_directory
+import shutil
+
+output_dir = 'output/'
+data_dir = './data/'
+make_directory(output_dir)
+shutil.copyfile(data_dir+'tanga2dxz.msh', output_dir+'mesh2dxz.msh')
+shutil.copyfile(data_dir+'tanga2dxz_show.msh', output_dir+'mesh2dxz_show.msh')
+
+
+domain = slim3d.Domain(data_dir+'tanga2dxz.msh', data_dir+'periodicMesh.txt', reference_density=997.22441884)
+
+equations = slim3d.Slim3d_equations(domain, temperature=True)
+equations.set_implicit_vertical_diffusion(True)
+#equations.set_implicit_vertical_advection(True)
+equations.set_vertical_viscosity('gotm')
+equations.set_vertical_diffusivity('gotm')
+equations.set_horizontal_viscosity('smagorinsky')
+equations.set_additional_artificial_horizontal_shear(data_dir+'additional_shear/additional_shear.idx')
+#equations.set_horizontal_diffusivity('okubo')
+#equations.set_lax_friedrichs(.1)
+equations.set_bottom_friction(True, z0B=0.005, z0S=0.02)
+equations.set_limiter(False)
+equations.set_initial_temperature('netcdf', (data_dir+'initial_temperature.nc','T_init'))
+equations.set_initial_salinity('vertical_gradient', None, 0, 0.)
+equations.set_wind_stress('stress', (data_dir+'windStress.nc','u'), (data_dir+'windStress.nc','v') )
+equations.set_coriolis((data_dir+'coriolis.nc', 'coriolis'))
+equations.set_boundary_coast('coast')
+equations.set_conservative_ale(False)#True)
+
+
+time_loop = slim3d.Loop(equations, 
+        time_step=600,#1800,
+        export_time_step=3600,#86400, 
+        final_time=86400*50, 
+        output_directory=output_dir)
+time_loop.export_elevation()
+time_loop.export_temperature()
+time_loop.export_z_coordinate()
+time_loop.export_uv(False)
+time_loop.export_uv2d(False)
+time_loop.export_w()
+time_loop.export_rho()
+
+#time_loop.loop()
+time_loop.setup()
+while time_loop.get_time() < time_loop.final_time:
+    time_loop.advance_one_time_step()
+    time_loop.check_sanity()
+    #print('******************')
+    #print(time_loop.equations._slimSolver.getUseConservativeALE())
+    if time_loop.check_export():
+        time_loop.print_iter_info()
+        [relVolErr, relSErr, relTErr] = time_loop.check_mass_conservation()
+        [absSDev, absTDev] = time_loop.check_tracer_consistency()
+        time_loop.export_fields()
+time_loop.terminate()

benchmarks/swe3d/tartinville/tartinville.py

-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-from dgpy import *
-from shutil import copyfile
-import math, os, inspect, time
-
-meshFile = 'square_20layers.msh'
-
-odir = 'output'
-if Msg.GetCommRank()==0 and not os.path.exists(odir) :
-  os.mkdir(odir)
-copyfile(meshFile,odir + '/mesh.msh')
-pythonFile = inspect.getfile(inspect.currentframe())
-copyfile(pythonFile,odir + '/' + pythonFile)
-
-s = slim3dSolver(meshFile,["bottom_Sea"],["top_Sea"])
-s.setSolveS(1)
-s.setSolveT(0)
-s.setSolveSImplicitVerticalDiffusion(0)
-s.setSolveUVImplicitVerticalDiffusion(0)
-s.setSolveTurbulence(0)
-s.setAdvectTurbulentEnergy(0)
-s.setComputeBottomFriction(0)
-s.setFlagUVLimiter(1)
-s.setUseOptimalLimiter(1)
-s.setUseModeSplit(1)
-s.exportDirectory = odir
-
-def sinit(s, xyz) :
-  for i in range(xyz.size1()) :
-    d = math.sqrt(xyz(i,0) * xyz(i,0) + xyz(i,1) * xyz(i,1))
-    S = 33.75 + 1.1 * math.pow(d/3000., 8)
-    si = S if xyz(i, 2) > -10 and d < 3e3 else 34.85
-    s.set(i, 0, si)
-
-slim3dParameters.setRho0(1025)
-f = s.functions
-d = s.getDofs(True)
-f.SInitFunc = functionPython(1, sinit, [function.getCoordinates()])
-f.coriolisFunc = functionConstant(1.15e-4)
-f.rhoFunc = linearStateEquation(d.SDof.getFunction(), 0.78, 33.75)
-
-# .-------------------------------.
-# |       create equations        |
-# '-------------------------------'
-e = s.createEquations()
-
-eta2d = d.etaDof2d.getFunction()
-uvAv2d = d.uvAvDof2d.getFunction()
-eta0 = functionConstant(0)
-
-horMomEq = e.horMomEq
-horMomEq.setLaxFriedrichsFactor(0.0)
-#horMomBndSides = horMomEq.newBoundaryWallLinear(0.5)
-#horMomBndSides = horMomEq.newBoundaryWall() # this is bad
-#horMomBndZero = horMomEq.new0FluxBoundary()
-horMomBndSymm = horMomEq.newSymmetryBoundary('')
-horMomBndWall = horMomEq.newBoundaryWall()
-toReplaceHorMomOpen = VectorFunctorConst([eta2d])
-replaceByHorMomOpen = VectorFunctorConst([eta0])
-horMomBndOpen = horMomEq.newOutsideValueBoundaryGeneric("", toReplaceHorMomOpen, replaceByHorMomOpen)
-horMomEq.addBoundaryCondition('top_Sea',horMomBndSymm)
-horMomEq.addBoundaryCondition('bottom_Sea',horMomBndWall) # must be wall for bath!
-horMomEq.addBoundaryCondition('Wall',horMomBndOpen)
-
-if s.getSolveUVImplicitVerticalDiffusion() :
-  vertMomUEq = e.vertMomUEq
-  vertMomUBndZero = vertMomUEq.new0FluxBoundary()
-  vertMomUBndBottom = vertMomUEq.newBoundaryBottom(f.bottomFrictionFunc, f.zBotDistFunc)
-  #vertMomUBndSymm = vertMomUEq.newSymmetryBoundary('')
-  vertMomUEq.addBoundaryCondition(['Wall','top_Sea'], vertMomUBndZero)
-  vertMomUEq.addBoundaryCondition('bottom_Sea',vertMomUBndBottom)
-
-wEq = e.wEq
-wBndZero = wEq.new0FluxBoundary()
-wBndSymm = wEq.newSymmetryBoundary('')
-wEq.addBoundaryCondition(['Wall','top_Sea'],wBndSymm) # must be symm!!
-wEq.addBoundaryCondition('bottom_Sea',wBndZero)
-
-SEq = e.SEq
-SEq.setLaxFriedrichsFactor(0.0)
-SBndZero = SEq.new0FluxBoundary()
-SBndOpen = SEq.newInFluxBoundary(f.SInitFunc)
-SEq.addBoundaryCondition('top_Sea',SBndZero)
-SEq.addBoundaryCondition('bottom_Sea',SBndZero)
-SEq.addBoundaryCondition('Wall',SBndOpen)
-
-eta2d = d.etaDof2d.getFunction()
-uvAv2d = d.uvAvDof2d.getFunction()
-
-eta2dEq = e.eta2dEq
-#etaBndZero = eta2dEq.new0FluxBoundary()
-toReplaceEtaOpen = VectorFunctorConst([function.getSolution()])
-replaceByEtaOpen = VectorFunctorConst([eta0])
-etaBndOpen = eta2dEq.newOutsideValueBoundaryGeneric("", toReplaceEtaOpen, replaceByEtaOpen)
-eta2dEq.addBoundaryCondition('Wall',etaBndOpen)
-
-uv2dEq = e.uv2dEq
-#uv2dBndWall = uv2dEq.newBoundaryWall()
-toReplaceUV2dOpen = VectorFunctorConst([eta2d])
-replaceByUV2dOpen = VectorFunctorConst([eta0])
-uv2dBndOpen = uv2dEq.newOutsideValueBoundaryGeneric("", toReplaceUV2dOpen, replaceByUV2dOpen)
-uv2dEq.addBoundaryCondition('Wall',uv2dBndOpen)
-
-# list of fields to export
-export = []
-export.append(dgIdxExporter(d.SDof, odir + "/salinity"))
-export.append(dgIdxExporter(d.uvDof, odir + "/uv", True))
-export.append(dgIdxExporter(d.etaDof2d, odir + "/eta"))
-
-# .-------------------------------.
-# |     setup time integrator     |
-# '-------------------------------'
-t = slim3dTimeIntegratorPC(s)
-s.timeIntegrator = t
-
-rk = dgERK(e.eta2dEq, None, DG_ERK_EULER)
-dtRK = rk.computeInvSpectralRadius(d.etaDof2d)
-CFL = 1
-dt = CFL*dtRK
-if Msg.GetCommRank() == 0 :
-  print('Runge-Kutta dt :' + str(dtRK) + ' CFL:' + str(CFL) + ' dt:' + str(dt))
-# try to estimate dt3d
-maxUV = 0.9 # max advective velocity
-dt3d = dtRK * math.sqrt(9.81*10) / maxUV
-M = int(0.72*dt3d/dt)
-M = 15
-if Msg.GetCommRank() == 0 :
-  print('Mode split ratio M :' + str(M))
-
-s.timeIntegrator.setTimeStep(dt)
-s.timeIntegrator.set3dTimeStepRatio(M)
-s.timeIntegrator.setExportInterval(3600)
-s.timeIntegrator.setEndTime(12*12*3600) # 6 days
-s.timeIntegrator.setExportAtEveryTimeStep(False)
-
-tic = time.time()
-while s.timeIntegrator.getTime() < s.timeIntegrator.getEndTime() :
-  s.timeIntegrator.advanceOneTimeStep()
-  s.timeIntegrator.checkSanity()
-  s.timeIntegrator.checkTracerConsistency()
-  if s.timeIntegrator.checkExport() :
-    Msg.Info("export %i" % s.timeIntegrator.getExportCount())
-    for e in export :
-      e.exportIdx(s.timeIntegrator.getExportCount(), s.timeIntegrator.getTime())
-  Msg.Info("%5i [ %s ] normuv %6.12e clock %.1fs" % (s.timeIntegrator.getIter(),time.strftime("%d %b %Y %H:%M:%S", time.gmtime(s.timeIntegrator.getTime())), d.uvDof.norm(), time.time() - tic))
-  if s.timeIntegrator.checkFullExport() :
-    s.timeIntegrator.exportAll()
-s.timeIntegrator.terminate(0)
+# Tartinville test-case
+import numpy as np
+### Generate Mesh ###
+from dgpy.scripts import Common
+Common.genMesh("square.geo", 2)
+print('mesh done')
+### Prepro ###
+import slimPre
+mesh_file = 'square.msh'
+nPart = slimPre.dgpy.Msg.GetCommSize()
+if nPart > 1:
+    if slimPre.dgpy.Msg.GetCommRank() == 0:
+        slimPre.dgpy.dgMeshPartition(mesh_file, nPart)
+    slimPre.dgpy.Msg.Barrier()
+    mesh_file = mesh_file[:-4] + '_' + str(nPart) + '.msh'
+mesh2d = slimPre.Mesh(mesh_file)
+slimPre.write_file('bath_2d.nc', region=None, time=None, data=[('bath',20.)])
+slimPre.netcdf_to_msh(mesh_file,'bath_2d.nc','bath','bath')
+slimPre.extrude(mesh_file, ('bath_2d.nc','bath'), nb_layers = 20, mesh_file_name_out='box.msh')
+mesh_file = 'box.msh'
+mesh = slimPre.Mesh(mesh_file)
+region_global = slimPre.Region(mesh)
+xyz = region_global.coordinates
+slimPre.write_file('coriolis.nc', region=None, time=None, data=[('coriolis',1.15e-4)])
+d = np.sqrt(xyz[:,0] * xyz[:,0] + xyz[:,1] * xyz[:,1])
+S = 33.75 + 1.1 * np.power(d/3000., 8)
+initS = np.where((xyz[:, 2] > -10) * (d < 3e3), S, 34.85)
+slimPre.write_file('initSalt.nc', region=region_global, time=None, data=[('salt',initS)])
+rho0 = 1025.
+print('prepro done')
+import slim3d
+domain = slim3d.Domain('box.msh', reference_density=rho0)
+equations = slim3d.Slim3d_equations(domain,salinity=True)
+equations.set_implicit_vertical_diffusion(False)
+equations.set_limiter(True)
+#equations.set_linear_solver_2d(False)
+equations.set_coriolis(('coriolis.nc','coriolis'))
+equations.set_bottom_friction(False)
+equations.set_initial_salinity("netcdf",('initSalt.nc', 'salt'))
+#equations.set_boundary_open(['Wall'], salinity=('initSalt.nc', 'salt'))
+equations.set_boundary_coast(['Wall'])
+equations.set_linear_density("salinity", 33.75, 0.78)
+time_loop = slim3d.Loop(equations, 
+        time_step=600.,
+        export_time_step=3600, 
+        final_time=6*86400,  
+        output_directory='output')
+time_loop.export_elevation()
+time_loop.export_uv(True)
+time_loop.export_rho()
+time_loop.export_uv2d(True)
+time_loop.export_salinity()
+#time_loop.set_timer()
+time_loop.loop()

dgpy/scripts/slim3d.py

@@ -250,6 +250,7 @@ class Loop:
         self._odir = output_directory
         self._export_uv       = False
         self._export_w        = False
+        self._export_z        = True
         self._export_S        = False
         self._export_T        = False
         self._export_uvAv2d   = False
@@ -273,6 +274,7 @@ class Loop:
         self._export_uv = True
         self._export_uv_vect = vect
     def export_w(self): self._export_w = True
+    def export_z_coordinate(self): self._export_z = True
     def export_salinity(self): self._export_S = True
     def export_temperature(self): self._export_T = True
     def export_uv2d(self, vect=True):
@@ -295,6 +297,29 @@ class Loop:
     def set_timer(self):
         self._timer = True
 
+    def export_fields(self):
+        timeIntegrator = self._timeIntegrator
+        if dgpy.Msg.GetCommRank()==0 : dgpy.Msg.Info("export %i" % timeIntegrator.getExportCount())
+        d = self._slimSolver.getDofs()
+        if self._export_nuh:
+            d.nuhDof3d.interpolate(self._slimSolver.functions.nuh)
+        if self._export_uvAv2d:
+            timeIntegrator.exporter_computeUVMean()
+        if self._export_w:
+            timeIntegrator.exporter_computeW()
+        if self._export_rho:
+            timeIntegrator.exporter_computeRho()
+        if self._export_eta:
+            timeIntegrator.exporter_computeEta()
+        if self._export_z:
+            timeIntegrator.exporter_computeZ()
+        for e in self._export :
+            e.exportIdx(timeIntegrator.getExportCount(), timeIntegrator.getTime() - self.initial_time)
+        if self._ncWriter3d:
+            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)
+
     def setup(self):
         slim3d_private.slim3d_setup(self)
         slimSolver = self._slimSolver
@@ -304,11 +329,12 @@ class Loop:
         dgpy.Msg.Barrier()
         d = slimSolver.getDofs()
         if self._export_uv:     self._export.append(dgpy.dgIdxExporter(d.uvDof, self._odir + "/uv", self._export_uv_vect))
-        if self._export_w:      self._export.append(dgpy.dgIdxExporter(d.wDof3d, self._odir + "/w"))
-        if self._export_rho:    self._export.append(dgpy.dgIdxExporter(d.rhoDof3d, self._odir + "/rho"))
+        if self._export_w:      self._export.append(dgpy.dgIdxExporter(d.wDof3dExport, self._odir + "/w"))
+        if self._export_z:      self._export.append(dgpy.dgIdxExporter(d.zDofExport, self._odir + "/z"))
+        if self._export_rho:    self._export.append(dgpy.dgIdxExporter(d.rhoDof3dExport, self._odir + "/rho"))
         if self._export_nuh:    self._export.append(dgpy.dgIdxExporter(d.nuhDof3d, self._odir + "/nuh"))
-        if self._export_uvAv2d: self._export.append(dgpy.dgIdxExporter(d.uvAvDof2d, self._odir + "/uvAv2d", self._export_uvAv2d_vect))
-        if self._export_eta:    self._export.append(dgpy.dgIdxExporter(d.etaDof2d, self._odir + "/eta"))
+        if self._export_uvAv2d: self._export.append(dgpy.dgIdxExporter(d.uvAvDof2dExport, self._odir + "/uvAv2d", self._export_uvAv2d_vect))
+        if self._export_eta:    self._export.append(dgpy.dgIdxExporter(d.etaDof2dExport, self._odir + "/eta"))
         if self._export_S: 
             if slimSolver.getSolveS():
                 self._export.append(dgpy.dgIdxExporter(d.SDof, self._odir + "/salinity"))
@@ -350,15 +376,7 @@ class Loop:
         timeIntegrator = self._timeIntegrator
         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())
-        if self._export_nuh:
-            d.nuhDof3d.interpolate(self._slimSolver.functions.nuh)
-        for e in self._export :
-            e.exportIdx(timeIntegrator.getExportCount(), timeIntegrator.getTime() - self.initial_time)
-        if self._ncWriter3d:
-            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)
+        self.export_fields()
         self._tic = slim_private.time.time()
 
     def get_time(self):
@@ -369,18 +387,6 @@ class Loop:
         self._timeIntegrator.checkSanity()
     def check_export(self):
         return self._timeIntegrator.checkExport()
-    def export_fields(self):
-        timeIntegrator = self._timeIntegrator
-        if dgpy.Msg.GetCommRank()==0 : dgpy.Msg.Info("export %i" % timeIntegrator.getExportCount())
-        d = self._slimSolver.getDofs()
-        if self._export_nuh:
-            d.nuhDof3d.interpolate(self._slimSolver.functions.nuh)
-        for e in self._export :
-            e.exportIdx(timeIntegrator.getExportCount(), timeIntegrator.getTime() - self.initial_time)
-        if self._ncWriter3d:
-            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)
     def check_full_export(self):
         return self._timeIntegrator.checkFullExport()
     def export_full(self):

dgpy/scripts/slim3d_private.py

@@ -206,7 +206,7 @@ def slim3d_setup(loop):
         eq._areaFunc = dgpy.functionPrecomputedExtrusion(slimSolver.extrusion(), 3, area2d.getFunction())
         if eq._additional_shear:
             eq._add_shear = slim_private._load_function(eq._additional_shear, slimSolver.groups3d)
-        f.nuh = dgpy.domeSmagorinsky(eq._hor_visc_fact, eq._hor_visc_max, d.uvGradDof.getFunction(), eq._areaFunc, eq._add_shear)
+        f.nuh = dgpy.domeSmagorinsky(eq._hor_visc_fact, eq._hor_visc_max, d.tmpuvGradDof.getFunction(), eq._areaFunc, eq._add_shear)
     
     if eq._horizontal_diffusivity == 'constant':
         f.kappahTotal = dgpy.functionConstant(eq._hor_diff_const)
@@ -222,14 +222,14 @@ def slim3d_setup(loop):
         f.nuvFunc = dgpy.functionConstant(eq._vertical_viscosity_constant_value)
     elif eq._vertical_viscosity == 'pacanowski_philander':
         paca_values = eq._vertical_viscosity_paca_values
-        f.nuvFunc = dgpy.nuPacanowskiPhilander(d.uvGradDof.getFunction(), d.rhoDof3d.getFunctionGradient(), paca_values[0], paca_values[1], paca_values[2], paca_values[3])
+        f.nuvFunc = dgpy.nuPacanowskiPhilander(d.tmpuvGradDof.getFunction(), d.tmprhoDof3d.getFunctionGradient(), paca_values[0], paca_values[1], paca_values[2], paca_values[3])
     if eq._vertical_diffusivity == 'constant':
         f.kappavFunc = dgpy.functionConstant(eq._vertical_diffusivity_constant_value)
     elif eq._vertical_diffusivity == 'pacanowski_philander':
         paca_values = eq._vertical_diffusivity_paca_values
         if not f.nuvFunc:
             dgpy.Msg.Fatal('To give a Pacanowski-Philander vertical diffusivity, a vertical viscosity is required')
-        f.kappavFunc = dgpy.kappaPacanowskiPhilander(d.uvGradDof.getFunction(), d.rhoDof3d.getFunctionGradient(), f.nuvFunc, paca_values[0], paca_values[1])
+        f.kappavFunc = dgpy.kappaPacanowskiPhilander(d.tmpuvGradDof.getFunction(), d.tmprhoDof3d.getFunctionGradient(), f.nuvFunc, paca_values[0], paca_values[1])
     
     if slimSolver.getComputeBottomFriction():
         f.z0BFunc = dgpy.functionConstant(eq._z0[0])
@@ -289,9 +289,10 @@ def slim3d_setup(loop):
         horMomEq.setBoundaryNeumann(openBnd.tag, openBnd.horizontal_momentum_flux_f)
 
     horMomBndWall = horMomEq.newBoundaryWall()
+    horMomBndWallVerticalBottom = horMomEq.newBoundaryWallVerticalBottom()
     horMomEq.addBoundaryCondition(eq._boundary_coast, horMomBndWall)
     horMomEq.addBoundaryCondition(bottomTags, horMomBndWall)
-    horMomEq.addBoundaryCondition('vertical_bottom', horMomBndWall)
+    horMomEq.addBoundaryCondition('vertical_bottom', horMomBndWallVerticalBottom)
     
     if slimSolver.getSolveUVImplicitVerticalDiffusion() :
         if f.windStressFunc :

dofContainer/dgDofContainer.cpp

@@ -780,23 +780,19 @@ void dgDofContainer::scale(double f)
 }
 
 void dgDofContainer::multiply(const dgDofContainer &x, const dgDofContainer &y){
-  fullMatrix<double> dataP, ghostDataP, dataX, ghostDataX, dataY, ghostDataY;
-  for (int iGroup = 0; iGroup < _groups->getNbElementGroups(); iGroup++) {
+  fullMatrix<double> dataP, dataX, dataY;
+  for (int iGroup = 0; iGroup < _groups->getNbElementGroups() + _groups->getNbGhostGroups(); iGroup++) {
     const dgGroupOfElements &group = *_groups->getElementGroup(iGroup);
     int nF = nField(group);
     if(nF != x.nField(group) || nF != y.nField(group))
       Msg::Fatal("vector size must match nbfields in dgDofContainer multiply(vector, vector)");
     for (size_t iElement = 0; iElement < group.getNbElements(); ++iElement) {
       getGroupProxy(iGroup).getBlockProxy(iElement, dataP);
-      getGroupProxy(iGroup).getBlockProxy(iElement, ghostDataP);
       x.getGroupProxy(iGroup).getBlockProxy(iElement, dataX);
-      x.getGroupProxy(iGroup).getBlockProxy(iElement, ghostDataX);
       y.getGroupProxy(iGroup).getBlockProxy(iElement, dataY);
-      y.getGroupProxy(iGroup).getBlockProxy(iElement, ghostDataY);
       for(int k=0; k < nF; k++) {
         for(int i=0; i < group.getNbNodes(); i++) {
           dataP(i, k) = dataX(i, k) * dataY(i, k);
-          ghostDataP(i, k) = ghostDataX(i, k) * ghostDataY(i, k);
         }
       }
     }
@@ -804,28 +800,23 @@ void dgDofContainer::multiply(const dgDofContainer &x, const dgDofContainer &y){
 }
 
 void dgDofContainer::divide(const dgDofContainer &x, const dgDofContainer &y, int iFieldY){
-  fullMatrix<double> dataP, ghostDataP, dataX, ghostDataX, dataY, ghostDataY;
-  for (int iGroup = 0; iGroup < _groups->getNbElementGroups(); iGroup++) {
+  fullMatrix<double> dataP, dataX, dataY;
+  for (int iGroup = 0; iGroup < _groups->getNbElementGroups() + _groups->getNbGhostGroups(); iGroup++) {
     const dgGroupOfElements &group = *_groups->getElementGroup(iGroup);
     int nF = nField(group);
     if( iFieldY == -1 && (nF != x.nField(group) || nF != y.nField(group)))
       Msg::Fatal("vector size must match nbfields in dgDofContainer divide(vector, vector)");
     for (size_t iElement = 0; iElement < group.getNbElements(); ++iElement) {
       getGroupProxy(iGroup).getBlockProxy(iElement, dataP);
-      getGroupProxy(iGroup).getBlockProxy(iElement, ghostDataP);
       x.getGroupProxy(iGroup).getBlockProxy(iElement, dataX);
-      x.getGroupProxy(iGroup).getBlockProxy(iElement, ghostDataX);
       y.getGroupProxy(iGroup).getBlockProxy(iElement, dataY);
-      y.getGroupProxy(iGroup).getBlockProxy(iElement, ghostDataY);
       for(int k=0; k < nF; k++) {
         for(int i=0; i < group.getNbNodes(); i++) {
           if (iFieldY == -1){
             dataP(i, k) = dataX(i, k) / dataY(i, k);
-            ghostDataP(i, k) = ghostDataX(i, k) / ghostDataY(i, k);
           }
           else{
             dataP(i, k) = dataX(i, k) / dataY(i, iFieldY);
-            ghostDataP(i, k) = ghostDataX(i, k) / ghostDataY(i, iFieldY);
           }
         }
       }

mesh/dgExtrusion.cpp

@@ -466,6 +466,19 @@ void dgExtrusion::_buildVertexMap()
       }
     }
   }
+  //reversemap
+  _vEdgeToDofDG.resize(getNbVerticals(false));
+  for (size_t i = 0; i < getNbVerticals(false); ++i) {
+    _vEdgeToDofDG[i].resize(getNbDGPointsInVertical(i));
+  }
+  for (size_t i = 0; i < _dofToVEdgeDG.size(); ++i) {
+    for (size_t j = 0; j < _dofToVEdgeDG[i].size(); ++j) {
+      for (size_t k = 0; k < _dofToVEdgeDG[i][j].size(); ++k) {
+        const std::pair<size_t, size_t> &v = _dofToVEdgeDG[i][j][k];
+        _vEdgeToDofDG[v.first][v.second] = {i, j, k};
+      }
+    }
+  }
 }
 
 

mesh/dgExtrusion.h

@@ -32,6 +32,8 @@ class dgExtrusion {
     std::vector<_nodeRef> _verticalTo2DNode;
     // maps each 2d node (iG2d, iEl2d, iN2d) to a vertical. Opposive of _verticalTo2DNode
     std::vector<std::vector<std::vector<size_t> > > _2dNodeToVertical;
+    // maps (iVertical, iLevel) to (iG,iE,iN) 3d
+    std::vector<std::vector<_nodeRef > > _vEdgeToDofDG;
     // generates unique int for FD nodes
     size_t _nbResidentCols;
     size_t _nbVerticals, _nbResidentVerticals;
@@ -129,6 +131,15 @@ class dgExtrusion {
     double innerRadius2d(size_t iGroup, size_t iElement) const;
     double area2d(size_t iGroup, size_t iElement) const;
     /** map 3D node to 2D node  (DG) **/
+    inline void map2d3d(size_t iGroup2, size_t iElement2, size_t iNode2, size_t iLevel, size_t &iGroup3, size_t &iElement3, size_t &iNode3) const
+    {
+      size_t iVertical = _2dNodeToVertical[iGroup2][iElement2][iNode2];
+      const _nodeRef &nr = _vEdgeToDofDG[iVertical][iLevel];
+      iElement3 = nr.elementId;
+      iGroup3 = nr.groupId;
+      iNode3 = nr.nodeId;
+    }
+    /** map 3D node to 2D node  (DG) **/
     inline void map3d2d(size_t iGroup3, size_t iElement3, size_t iNode3, size_t &iGroup2, size_t &iElement2, size_t &iNode2) const
     {
       size_t iVertical, iLevel;

modules/shallowWater/dgConservationLawShallowWater2d.cpp

@@ -136,10 +136,12 @@ void dgConservationLawShallowWater2d::gradPsiTerm(functorCache &cm, fullMatrix<d
       double u = sol(i,1);
       double v = sol(i,2);
       double A = movingBathFactor(i,0);//equals 1 if no WD
-      if (!_linear)
+      if (!_linear){
         H+=eta;
-      else if (_linearSolverFrom3D)
+      }
+      else if (_linearSolverFrom3D){
         H+=(*etaOld)(i,0);
+      }
       // flux_x
       val(i,0) = H*u /A;
       val(i,1) = g(i,0)*(1+rhoSurf(i,0)/_density)*eta + atmPress(i,0)/_density;
@@ -206,7 +208,7 @@ void dgConservationLawShallowWater2d::source(functorCache &cm, fullMatrix<double
         val(i,1) += (*nu)(i, 0)/H*(dHdx*(*solGradient)(i, 3)+dHdy*(*solGradient)(i, 4));
         val(i,2) += (*nu)(i, 0)/H*(dHdx*(*solGradient)(i, 6)+dHdy*(*solGradient)(i, 7));
       }
-      if(!_linear){
+      if(!_linear && !_from3D){
         double div = (*solGradient)(i,3) + (*solGradient)(i,7);
         double bottomFriction = hypot(u,v)*quadraticDissipation(i,0);
         val(i,1) += u*(- bottomFriction + div);
@@ -326,15 +328,16 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
       double c, utM = 0.;
       if (_linear) {
         c = sqrt(g(i,0) * H);
-        if (_linearSolverFrom3D)
+        if (_linearSolverFrom3D){
           H += ((*etaOldL)(i,0) + (*etaOldR)(i,0))/2;
+        }
       }
       else {
         H += etaM;
         vR = vR * hR/hMin;
         vL = vL * hL/hMin;
         utM = (vR + vL) /2;
-        c = sqrt(g(i,0) * H + fabs(unM));
+        c = sqrt(g(i,0) * H/* + fabs(unM)*/);
       }
       Fun  = - g(i,0) * (1+rhoSurf/_density) * etaM - pa/_density + c * (uR - uL) / 2;
       Feta = - H * unM + c * (etaR - etaL) / 2;
@@ -344,6 +347,13 @@ void dgConservationLawShallowWater2d::riemann(functorCache &cache,fullMatrix<dou
         //Feta /= Amb; TODO verify if Wetting and drying and Lax Friedrichs !
         Fut += - unM * utM + c * (vR - vL) / 2; // jump only if advection
       }
+      /*
+      if (cache.interfaceGroup()->nConnection() == 2) {
+        int el0 = cache.interfaceGroup()->element(cache.interfaceId(i/cache.nPointByElement()),0).num();
+        int el1 = cache.interfaceGroup()->element(cache.interfaceId(i/cache.nPointByElement()),1).num();
+        if (el0 == 560 && el1 == 578)
+          printf("flux * H in 2d is %.16g eta is %.16g %.16g H (%.16g %.16g) is %.16g u is %.16g %.16g (%.16g)\n", Feta, solL(i,0), solR(i,0), H, hL, hR, uL/hMin*hL, uR/hMin*hR, uR);
+      }*/
     }
     else {
       if (_from3D)

modules/slim3d/dgConservationLawSW3dContinuity.cpp

@@ -36,34 +36,26 @@ void dgConservationLawSW3dContinuity::riemann(functorCache &cache, fullMatrix<do
   const fullMatrix<double> &etaOut = cache.get(*_eta, 1);
   const fullMatrix<double> &UIn = cache.get(*_uv, 0);