Merge pull request #20 from ulgltas/adrien

CUPyDO v1.2.3 (flow compatibility update)

cupydo/interfaces/Flow.py

@@ -58,10 +58,19 @@ def __initFlow(self, p, _nthreads):
         import flow
         import tbox
         import tbox.gmsh as gmsh
+        from tbox.solvers import LinearSolver
 
         # basic checks
         if p['Dim'] != 2 and p['Dim'] != 3:
-            raise Exception('Problem dimension should be 2 or 3, but ' + p['Dim'] + ' was given!\n')
+            raise RuntimeError('Problem dimension should be 2 or 3, but ' + p['Dim'] + ' was given!\n')
+        if p['Dim'] == 2:
+            if 'AoS' in p and p['AoS'] != 0:
+                raise RuntimeError('Angle of sideslip (AoS) should be zero for 2D problems!\n')
+            if 'Symmetry' in p:
+                raise RuntimeError('Symmetry boundary condition cannot be used for 2D problems!\n')
+        else:
+            if 'AoS' in p and p['AoS'] != 0 and 'Symmetry' in p:
+                raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n')
         # basic config
         if p['Format'] == 'vtk':
             try:   
@@ -86,10 +95,11 @@ def __initFlow(self, p, _nthreads):
             for i in range(0, len(p['Wakes'])):
                 mshCrck = tbox.MshCrack(self.msh, p['Dim'])
                 mshCrck.setCrack(p['Wakes'][i])
+                mshCrck.addBoundaries([p['Fluid'], p['Farfield'][-1], p['Wings'][i]])
                 if 'Fuselage' in p:
-                    mshCrck.addBoundaries([p['Fluid'], p['Symmetry'], p['Farfield'][-1], p['Fuselage'], p['Wings'][i]])
-                else:
-                    mshCrck.addBoundaries([p['Fluid'], p['Symmetry'], p['Farfield'][-1], p['Wings'][i]])
+                    mshCrck.addBoundaries([p['Fuselage']])
+                if 'Symmetry' in p:
+                    mshCrck.addBoundaries([p['Symmetry']])
                 mshCrck.setExcluded(p['WakeTips'][i])
                 mshCrck.run()
         tbox.GmshExport(self.msh).save(self.msh.name)
@@ -118,13 +128,14 @@ def __initFlow(self, p, _nthreads):
 
         # initialize the problem
         p['AoA'] = p['AoA']*np.pi/180 # convert to radians
-        phiInfFun = flow.F0PsPhiInf(p['Dim'], p['AoA'])
-        if p['Dim'] == 2:
-            velInfFun = tbox.Fct1C(-np.cos(p['AoA']), -np.sin(p['AoA']), 0.)
+        if 'AoS' in p:
+            p['AoS'] = p['AoS']*np.pi/180
         else:
-            velInfFun = tbox.Fct1C(-np.cos(p['AoA']), 0., -np.sin(p['AoA']))
+            p['AoS'] = 0.
+        phiInfFun = flow.F0PsPhiInf(p['Dim'], p['AoA'], p['AoS'])
+        velInfFun = flow.F1ElVi(p['Dim'], p['AoA'], p['AoS'])
         self.dynP = p['P_dyn']    
-        pbl = flow.Problem(self.msh, p['Dim'], p['AoA'], p['M_inf'], p['S_ref'], p['c_ref'], p['x_ref'], p['y_ref'], p['z_ref'])
+        pbl = flow.Problem(self.msh, p['Dim'], p['AoA'], p['AoS'], p['M_inf'], p['S_ref'], p['c_ref'], p['x_ref'], p['y_ref'], p['z_ref'])
 
         # add medium
         if p['M_inf'] == 0:
@@ -157,21 +168,36 @@ def __initFlow(self, p, _nthreads):
             for i in range(0, len(p['Wakes'])):
                 pbl.add(flow.Wake(self.msh, [p['Wakes'][i], p['Wakes'][i]+'_', p['Fluid'], p['TeTips'][i]]))
 
-        # initialize the solver
+        # initialize the linear (inner) solver
+        if p['LSolver'] == 'Pardiso':
+            linsol = LinearSolver().pardiso()
+        elif p['LSolver'] == 'GMRES':
+            linsol = tbox.Gmres()
+            linsol.setFillFactor(p['G_fill'])
+            linsol.setRestart(p['G_restart'])
+            linsol.setTolerance(p['G_tol'])
+        elif p['LSolver'] == 'MUMPS':
+            linsol = LinearSolver().mumps()
+        elif p['LSolver'] == 'SparseLU':
+            linsol = tbox.SparseLu()
+        else:
+            raise RuntimeError('Available linear solvers: Pardiso, GMRES, MUMPS or SparseLU, but ' + p['LSolver'] + ' was given!\n')
+        # initialize the nonlinear (outer) solver
         if p['NSolver'] == 'Picard':
-            self.solver = flow.Picard(pbl)
+            self.solver = flow.Picard(pbl, linsol)
             self.solver.relax = p['Relaxation']
         elif p['NSolver'] == 'Newton':
-            self.solver = flow.Newton(pbl)
+            self.solver = flow.Newton(pbl, linsol)
             self.solver.lsTol = p['LS_tol']
             self.solver.maxLsIt = p['Max_it_LS']
             self.solver.avThrsh = p['AV_thrsh']
         else:
-            raise RuntimeError('Available nonlinear solver type: Picard or Newton, but ' + p['NSolver'] + ' was given!\n')
+            raise RuntimeError('Available nonlinear solvers: Picard or Newton, but ' + p['NSolver'] + ' was given!\n')
         self.solver.nthreads = _nthreads
         self.solver.relTol = p['Rel_tol']
         self.solver.absTol = p['Abs_tol']
         self.solver.maxIt = p['Max_it']
+        print 'Linear solver: ', linsol
         print "Number of threads: ", self.solver.nthreads
         print "Maximum number of iterations: ", self.solver.maxIt
         print "Objective relative residual: ", self.solver.relTol
@@ -207,9 +233,9 @@ def getNodalInitialPositions(self):
         z0 = np.zeros(self.nPhysicalNodes)
         for i in range(self.boundary.nodes.size()):
             n = self.boundary.nodes[i]               
-            x0[i] = n.pos.x[0]
-            y0[i] = n.pos.x[1]
-            z0[i] = n.pos.x[2]
+            x0[i] = n.pos[0]
+            y0[i] = n.pos[1]
+            z0[i] = n.pos[2]
 
         return (x0, y0, z0)
 
@@ -226,9 +252,9 @@ def applyNodalDisplacements(self, dx, dy, dz, dx_nM1, dy_nM1, dz_nM1, haloNodesD
         """
         self.mshDef.savePos()
         for i in range(self.boundary.nodes.size()):
-            self.boundary.nodes[i].pos.x[0] = self.nodalInitPosX[i] + dx[i]
-            self.boundary.nodes[i].pos.x[1] = self.nodalInitPosY[i] + dy[i]
-            self.boundary.nodes[i].pos.x[2] = self.nodalInitPosZ[i] + dz[i]
+            self.boundary.nodes[i].pos[0] = self.nodalInitPosX[i] + dx[i]
+            self.boundary.nodes[i].pos[1] = self.nodalInitPosY[i] + dy[i]
+            self.boundary.nodes[i].pos[2] = self.nodalInitPosZ[i] + dz[i]
 
     def meshUpdate(self, nt):
         """Deform the mesh using linear elasticity equations

tests/Flow_Metafor/staticAgard_fluid.py

@@ -40,6 +40,7 @@ def getParams():
     p['y_ref'] = 0. # Reference point for moment computation (y)
     p['z_ref'] = 0. # Reference point for moment computation (z)
     # Numerical
+    p['LSolver'] = 'Pardiso' # Linear solver (Pardiso, GMRES, MUMPS or SparseLU)
     p['NSolver'] = 'Newton' # Noninear solver type (Picard or Newton)
     p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual
     p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual

tests/Flow_Metafor/staticDiamond_fluid.py

@@ -31,6 +31,7 @@ def getParams():
     p['y_ref'] = 0. # Reference point for moment computation (y)
     p['z_ref'] = 0. # Reference point for moment computation (z)
     # Numerical
+    p['LSolver'] = 'Pardiso' # Linear solver (Pardiso, GMRES, MUMPS or SparseLU)
     p['NSolver'] = 'Newton' # Noninear solver type (Picard or Newton)
     p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual
     p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual

tests/Flow_Modal/staticAgard_fluid.py

@@ -40,6 +40,10 @@ def getParams():
     p['y_ref'] = 0. # Reference point for moment computation (y)
     p['z_ref'] = 0. # Reference point for moment computation (z)
     # Numerical
+    p['LSolver'] = 'GMRES' # Linear solver (Pardiso, GMRES, MUMPS or SparseLU)
+    p['G_fill'] = 2 # Fill-in factor for GMRES preconditioner
+    p['G_tol'] = 1e-5 # Tolerance for GMRES
+    p['G_restart'] = 50 # Restart for GMRES
     p['NSolver'] = 'Newton' # Noninear solver type (Picard or Newton)
     p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual
     p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual

tests/Flow_RBM/staticNaca_fluid.py

@@ -33,6 +33,7 @@ def getParams():
     p['y_ref'] = 0. # Reference point for moment computation (y)
     p['z_ref'] = 0. # Reference point for moment computation (z)
     # Numerical
+    p['LSolver'] = 'Pardiso' # Linear solver (Pardiso, GMRES, MUMPS or SparseLU)
     p['NSolver'] = 'Newton' # Noninear solver type (Picard or Newton)
     p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual
     p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual