added substepping for custom integration scheme

.travis.yml

@@ -29,7 +29,7 @@ before_install:
 install:
 - git lfs fetch
 - git lfs checkout
-- git clone https://github.com/erdc-cm/proteus
+- git clone https://github.com/erdc/proteus
 - cd proteus
 - make hashdist
 - make stack

2d/floatingStructures/floatingCaissonAddedMass/MeshRefinement.py

@@ -0,0 +1,75 @@
+def geometry_to_gmsh(domain):
+    import py2gmsh
+    from py2gmsh.Mesh import *
+    from py2gmsh.Entity import *
+    from py2gmsh.Fields import *
+    self = domain
+    lines_dict = {}
+
+    mesh = Mesh()
+
+    if self.boundaryTags:
+        for tag, flag in self.boundaryTags.items():
+            phys = PhysicalGroup(nb=flag, name=tag)
+            mesh.addGroup(phys)
+
+    for i, v in enumerate(self.vertices):
+        if domain.nd == 2:
+            p = Point([v[0], v[1], 0.])
+        else:
+            p = Point(v)
+        mesh.addEntity(p)
+        g = mesh.groups.get(self.vertexFlags[i])
+        if g:
+            g.addEntity(p)
+    nb_points = i+1
+    for i in range(nb_points):
+        lines_dict[i] = {}
+
+    for i, s in enumerate(self.segments):
+        lines_dict[s[0]][s[1]] = i
+        l = Line([mesh.points[s[0]+1], mesh.points[s[1]+1]])
+        mesh.addEntity(l)
+        g = mesh.groups.get(self.segmentFlags[i])
+        if g:
+            g.addEntity(l)
+
+    for i, f in enumerate(self.facets):
+        if self.nd == 3 or (self.nd == 2 and i not in self.holes_ind):
+            lineloops = []
+            for j, subf in enumerate(f):
+                lineloop = []
+                # vertices in facet
+                for k, ver in enumerate(subf):
+                    if ver in lines_dict[subf[k-1]].keys():
+                        lineloop += [lines_dict[subf[k-1]][ver]+1]
+                    elif subf[k-1] in lines_dict[ver].keys():
+                        # reversed
+                        lineloop += [(lines_dict[ver][subf[k-1]]+1)]
+                    else:
+                        l = Line([mesh.points[subf[k-1]+1], mesh.points[ver+1]])
+                        mesh.addEntity(l)
+                        lineloop += [l.nb]
+                ll = LineLoop(mesh.getLinesFromIndex(lineloop))
+                mesh.addEntity(ll)
+                lineloops += [ll.nb]
+            s = PlaneSurface([mesh.lineloops[loop] for loop in lineloops])
+            mesh.addEntity(s)
+            g = mesh.groups.get(self.facetFlags[i])
+            if g:
+                g.addEntity(s)
+
+    for i, V in enumerate(self.volumes):
+        surface_loops = []
+        hole_loops = []
+        for j, sV in enumerate(V):
+            sl = SurfaceLoop(mesh.getSurfacesFromIndex((np.array(sV)+1).tolist()))
+            mesh.addEntity(sl)
+            surface_loops += [sl]
+        vol = Volume(surface_loops)
+        mesh.addEntity(vol)
+        g = mesh.groups.get(self.regionFlags[i])
+        if g:
+            g.addEntity(vol)
+
+    return mesh

2d/floatingStructures/floatingCaissonAddedMass/added_massIB_n.py

@@ -0,0 +1,56 @@
+from proteus.default_n import *
+import added_massIB_p as physics
+from proteus import (StepControl,
+                     TimeIntegration,
+                     NonlinearSolvers,
+                     LinearSolvers,
+                     LinearAlgebraTools)
+from proteus.mprans import AddedMass
+from proteus import Context
+
+ct = Context.get()
+domain = ct.domain
+nd = ct.domain.nd
+mesh = domain.MeshOptions
+
+elementQuadrature = ct.elementQuadrature
+elementBoundaryQuadrature = ct.elementBoundaryQuadrature
+
+triangleOptions = mesh.triangleOptions
+
+femSpaces = {0:ct.basis}
+
+stepController=StepControl.FixedStep
+
+numericalFluxType =  AddedMass.NumericalFlux
+
+matrix = LinearAlgebraTools.SparseMatrix
+
+
+multilevelLinearSolver = LinearSolvers.KSP_petsc4py
+levelLinearSolver      = LinearSolvers.KSP_petsc4py
+parallelPartitioningType = mesh.parallelPartitioningType
+nLayersOfOverlapForParallel = mesh.nLayersOfOverlapForParallel
+nonlinearSmoother = NonlinearSolvers.AddedMassNewton
+linearSmoother    = LinearSolvers.NavierStokesPressureCorrection
+
+linear_solver_options_prefix = 'am_'
+
+multilevelNonlinearSolver = NonlinearSolvers.AddedMassNewton
+levelNonlinearSolver = NonlinearSolvers.AddedMassNewton
+
+#linear solve rtolerance
+
+linTolFac = 0.0
+l_atol_res = 1.0e-10#0.001*ct.am_nl_atol_res
+tolFac = 0.0
+nl_atol_res = ct.am_nl_atol_res
+
+nonlinearSolverConvergenceTest = 'rits'
+levelNonlinearSolverConvergenceTest = 'rits'
+linearSolverConvergenceTest             = 'r-true'
+maxNonlinearIts =1
+maxLineSearches=0
+periodicDirichletConditions=None
+conservativeFlux=None
+auxiliaryVariables=[ct.system.ProtChAddedMass]

2d/floatingStructures/floatingCaissonAddedMass/added_massIB_p.py

@@ -0,0 +1,44 @@
+from proteus.default_p import *
+from proteus.mprans import AddedMass
+from proteus import Context
+name = "added_mass"
+ct = Context.get()
+domain = ct.domain
+nd = domain.nd
+mesh = domain.MeshOptions
+genMesh = mesh.genMesh
+movingDomain = ct.movingDomain
+T = ct.T
+
+LevelModelType = AddedMass.LevelModel
+
+coefficients = AddedMass.Coefficients(nd=nd,
+                                      V_model=4,
+                                      barycenters=domain.barycenters,
+                                      flags_rigidbody=ct.flags_rigidbody)
+
+def getDBC_am(x,flag):
+#    if flag == ct.domain.boundaryTags['tank2D1_y+']:
+#        return lambda x,t: 0.0
+#    else:
+    return None
+
+dirichletConditions = {0:getDBC_am}
+
+advectiveFluxBoundaryConditions = {}
+
+def getFlux_am(x, flag):
+    #the unit rigid motions will be applied internally
+    #leave this set to zero
+#    if flag == ct.domain.boundaryTags['tank2D1_y+']:
+#        return None
+#    else:
+    return lambda x,t: 0.0
+
+diffusiveFluxBoundaryConditions = {0: {0:getFlux_am}}
+
+class dp_IC:
+    def uOfXT(self, x, t):
+        return 0.0
+
+initialConditions = {0: dp_IC()}

2d/floatingStructures/floatingCaissonAddedMass/added_mass_n.py

@@ -42,7 +42,7 @@
 #linear solve rtolerance
 
 linTolFac = 0.0
-l_atol_res = ct.am_nl_atol_res
+l_atol_res = 1.0e-10#0.001*ct.am_nl_atol_res
 tolFac = 0.0
 nl_atol_res = ct.am_nl_atol_res
 

2d/floatingStructures/floatingCaissonAddedMass/added_mass_p.py

@@ -18,15 +18,21 @@
                                       flags_rigidbody=ct.flags_rigidbody)
 
 def getDBC_am(x,flag):
+#    if flag == ct.domain.boundaryTags['tank2D1_y+']:
+#        return lambda x,t: 0.0
+#    else:
     return None
 
 dirichletConditions = {0:getDBC_am}
 
 advectiveFluxBoundaryConditions = {}
 
 def getFlux_am(x, flag):
-    #the unit rigid motions will applied internally
+    #the unit rigid motions will be applied internally
     #leave this set to zero
+#    if flag == ct.domain.boundaryTags['tank2D1_y+']:
+#        return None
+#    else:
     return lambda x,t: 0.0
 
 diffusiveFluxBoundaryConditions = {0: {0:getFlux_am}}