Merge pull request #29 from jefalon/dev

Fixed mpi optimizations

.github/workflows/CI_every_PR.yml

@@ -32,11 +32,23 @@ jobs:
           bash ./install.sh test-environment
           conda activate test-environment
         
-      - name: Run all tests within WindSE
+      - name: Run all documented demos within WindSE
         shell: pwsh
         run: |
           conda activate test-environment
-          pytest -v --cov=windse tests
+          pytest -sv --cov=windse tests/test_demos.py
+
+      - name: Run all regression tests within WindSE
+        shell: pwsh
+        run: |
+          conda activate test-environment
+          pytest -sv --cov=windse --cov-append tests/test_regression.py
+        
+      - name: Run parallel regression tests within WindSE
+        shell: pwsh
+        run: |
+          conda activate test-environment
+          pytest -sv --cov=windse --cov-append tests/test_regression_parallel.py
       
       # Run coveralls
       - name: Run coveralls

install.sh

@@ -8,7 +8,7 @@ conda create -y --name $1
 conda activate $1
 
 ### Install conda-forge dependencies 
-conda install -y -c conda-forge fenics=2019.1.0=py38_9 dolfin-adjoint matplotlib scipy=1.4.1 slepc mshr hdf5 pyyaml memory_profiler pytest pytest-cov coveralls
+conda install -y -c conda-forge fenics=2019.1.0=py38_9 dolfin-adjoint matplotlib scipy=1.4.1 slepc mshr hdf5 pyyaml memory_profiler pytest pytest-cov pytest-mpi coveralls
 
 ### Install the new tsfc compiler
 pip install git+https://github.com/blechta/[email protected]

tests/test_regression.py

@@ -12,7 +12,7 @@
 home_path = os.getcwd()
 reg_path = "../9-Regression/"
 
-### Get Yaml Files ###
+### Get Yaml Files to be run in serial (all of them) ###
 yaml_files = sorted(pathlib.Path(__file__, reg_path).resolve().glob('*.yaml'))
 
 ### Import the tolerances ###
@@ -29,7 +29,7 @@
 ### Run Demo Yaml Files
 @pytest.mark.parametrize('yaml_file', yaml_files, ids=lambda yaml_file: yaml_file.parts[-2]+"/"+yaml_file.parts[-1])
 def test_yaml_execution(yaml_file):
-
+    
     ### Filter out some benign numpy warnings ###
     warnings.filterwarnings("ignore", message="numpy.dtype size changed")
     warnings.filterwarnings("ignore", message="numpy.ufunc size changed")
@@ -39,8 +39,9 @@ def test_yaml_execution(yaml_file):
     os.chdir(folder)
     from windse_driver import driver
     driver.run_action(params_loc=yaml_file.as_posix())
+    # comm.Barrier()
     os.chdir(home_path)
-
+    
     ### Grab the name of the run ###
     _, yaml_name = os.path.split(yaml_file)
     yaml_name = yaml_name.split(".")[0]

windse/ObjectiveFunctions.py

@@ -35,8 +35,8 @@ def CalculatePowerFunctional(solver,inflow_angle = 0.0):
                 tf1 = solver.problem.farm.actuator_disks_list[i] * cos(yaw)**2
                 tf2 = solver.problem.farm.actuator_disks_list[i] * sin(yaw)**2
                 tf3 = solver.problem.farm.actuator_disks_list[i] * 2.0 * cos(yaw) * sin(yaw)
-                tf = tf1*solver.u_k[0]**2+tf2*solver.u_k[1]**2+tf3*solver.u_k[0]*solver.u_k[1]
-                J_list[i+1] = assemble(dot(tf,solver.u_k)*dx,**solver.extra_kwarg)
+                tf = tf1*solver.problem.u_k[0]**2+tf2*solver.problem.u_k[1]**2+tf3*solver.problem.u_k[0]*solver.problem.u_k[1]
+                J_list[i+1] = assemble(dot(tf,solver.problem.u_k)*dx,**solver.extra_kwarg)
             else:
                 print("WARNING: missing individual turbine actuator disk, only able to report full farm power")
 

windse/OptimizationManager.py

@@ -252,7 +252,7 @@ def DebugOutput(self):
             ### Output gradient ###
             if hasattr(self,"gradient"):
                 for i, d in enumerate(self.gradients):
-                    self.tag_output("grad_"+self.names[i],float(d), collective_output='sum')
+                    self.tag_output("grad_"+self.names[i],float(d))
 
             ### TODO: Output taylor convergence data
             if hasattr(self,"conv_rate"):
@@ -423,10 +423,10 @@ def Gradient(self):
             ctl_val = float(self.controls[i].values())
             # d_out = str(self.names[i] + "    " +repr(ctl_val)+ "    " +repr(float(d)))
 
-            d_format = np.float64(d)
-            self.params.comm.Gather(d_format, d_global, root=0)
-            d_sum = np.sum(d_global)
-            d_out = '%12s: %12.5e, %22.15e' % (self.names[i], ctl_val, d_sum)
+            # d_format = np.float64(d)
+            # self.params.comm.Gather(d_format, d_global, root=0)
+            # d_sum = np.sum(d_global)
+            d_out = '%12s: %12.5e, %22.15e' % (self.names[i], ctl_val, d)
 
             # print('Rank %d, %s' % (self.params.rank, d_out))
             self.fprint(d_out)

windse/ParameterManager.py

@@ -206,6 +206,10 @@ def Load(self, loc,updated_parameters=[]):
             for sub in subfolder_list:
                 os.makedirs('%s/%s' % (self.folder, sub), exist_ok=True)
 
+            os.makedirs('%s/data/alm' % (self.folder), exist_ok=True)
+            os.makedirs('%s/data/alm/rotor_force' % (self.folder), exist_ok=True)
+            os.makedirs('%s/data/alm/angle_of_attack' % (self.folder), exist_ok=True)
+
         # Wait until rank 0 has created the directory structure
         self.comm.barrier()
 

windse/ProblemManager.py

@@ -155,8 +155,10 @@ def ComputeTurbineForce(self,u,inflow_angle,simTime=0.0):
             ### create output files for alm data ###
             force_folder = self.params.folder+"data/alm/rotor_force/"
             aoa_folder = self.params.folder+"data/alm/angle_of_attack/"
-            if not os.path.exists(aoa_folder): os.makedirs(aoa_folder)
-            if not os.path.exists(force_folder): os.makedirs(force_folder)
+            # if not os.path.exists(aoa_folder): os.makedirs(aoa_folder)
+            # if not os.path.exists(force_folder): os.makedirs(force_folder)
+            if not os.path.exists(aoa_folder) and self.params.rank==0: os.makedirs(aoa_folder)
+            if not os.path.exists(force_folder) and self.params.rank==0: os.makedirs(force_folder)
             self.aoa_files = []
             self.force_files = []
             for i in range(self.farm.numturbs):

windse/SolverManager.py

@@ -815,14 +815,17 @@ def Solve(self,iter_val=0):
 
         # solver_1 = PETScKrylovSolver('gmres', 'jacobi')
         solver_1 = PETScKrylovSolver('gmres', 'default')
+        # solver_1 = PETScKrylovSolver('default', 'default')
         solver_1.set_operator(A1)
 
         # solver_2 = PETScKrylovSolver('gmres', 'petsc_amg')
         solver_2 = PETScKrylovSolver('gmres', 'hypre_amg')
+        # solver_2 = PETScKrylovSolver('default', 'default')
         solver_2.set_operator(A2)
 
         # solver_3 = PETScKrylovSolver('cg', 'jacobi')
         solver_3 = PETScKrylovSolver('gmres', 'default')
+        # solver_3 = PETScKrylovSolver('default', 'default')
         solver_3.set_operator(A3)
 
         pr = cProfile.Profile()
@@ -1089,6 +1092,7 @@ def save_adj(adj):
     def SaveTimeSeries(self, simTime, simIter=None):
 
         self.DebugOutput(simTime,simIter)
+        ### TODO THIS NEED TO BE CLEAN TO ACCOUNT FOR DISKS
 
         if hasattr(self.problem,"tf_save"):
             self.problem.tf_save.vector()[:] = 0

windse/WindFarmManager.py

@@ -1084,7 +1084,7 @@ def init_constant_alm_terms(problem):
                 problem.min_dist.append(dist)
 
             # Create a Constant "wrapper" to enable dolfin to track mpi_u_fluid
-            problem.mpi_u_fluid_constant = Constant(np.zeros((problem.farm.numturbs, 3*3*problem.num_blade_segments)))
+            problem.mpi_u_fluid_constant = Constant(np.zeros((problem.farm.numturbs, 3*3*problem.num_blade_segments)),name="mpi_u_fluid")
 
 
         def init_unsteady_alm_terms(problem):
@@ -1162,7 +1162,7 @@ def init_mpi_alm(problem):
                     # If this rank doesn't own that point, an error will occur,
                     # in which case zeros should be reported
                     try:
-                        fn_val = problem.u_k1(xi, yi, zi)
+                        fn_val = problem.u_k1(np.array([xi, yi, zi]))
                         mpi_u_fluid[k, 3*j:3*j+3] = fn_val
                         mpi_u_fluid_count[k, 3*j:3*j+3] = [1, 1, 1]
                     except:
@@ -1231,7 +1231,7 @@ def finalize_mpi_alm(problem):
         mpi_u_fluid = init_mpi_alm(problem)
 
         # Populate the Constant "wrapper" with the velocity values to enable dolfin to track mpi_u_fluid
-        problem.mpi_u_fluid_constant.assign(Constant(mpi_u_fluid))
+        problem.mpi_u_fluid_constant.assign(Constant(mpi_u_fluid,name="temp_u_f"))
 
         # Call the ALM function for each turbine individually
         alm_output_list = []

windse/dolfin_adjoint_helper.py

@@ -91,8 +91,9 @@ def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx, prepar
             adj = self.ground(x,y,dy=1)
             # adj = (self.ground(x,y+h)-self.ground(x,y-h))/(2*h)
 
+        adj_output = adj_input * adj
 
-        return adj_input * adj
+        return adj_output
 
 
 def RadialChordForce(r,chord):
@@ -268,6 +269,12 @@ def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx, prepar
         for i in range(3):
             adj_output += np.inner(adj_inputs[i].get_local(self.all_ids),prepared[name][i][:,turb_idx])
 
+        adj_output = np.float64(adj_output)
+        recv_buff = np.zeros(1, dtype=np.float64)
+        self.farm.params.comm.Allreduce(adj_output, recv_buff)
+        adj_output = recv_buff
+
+        # adj_output = dolfin.MPI.sum(dolfin.MPI.comm_world,adj_output)
         return np.array(adj_output)
 
 
@@ -912,12 +919,24 @@ def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx, prepar
         elif name == "yaw":
             comp_vec = prepared[name]
             adj_vec = adj_inputs[0].get_local()
-            adj_output = np.inner(comp_vec, adj_vec)
+            adj_output = np.float64(np.inner(comp_vec, adj_vec))
+
+            recv_buff = np.zeros(1, dtype=np.float64)
+            self.problem.params.comm.Allreduce(adj_output, recv_buff)
+            adj_output = recv_buff
+
+            # adj_output = dolfin.MPI.sum(dolfin.MPI.comm_world,adj_output)
             return np.array(adj_output)
         else:
             comp_vec = prepared[name][:, segment_index]
             adj_vec = adj_inputs[0].get_local()
-            adj_output = np.inner(comp_vec, adj_vec)
+            adj_output = np.float64(np.inner(comp_vec, adj_vec))
+
+            recv_buff = np.zeros(1, dtype=np.float64)
+            self.problem.params.comm.Allreduce(adj_output, recv_buff)
+            adj_output = recv_buff
+
+            # adj_output = dolfin.MPI.sum(dolfin.MPI.comm_world,adj_output)
             return np.array(adj_output)
 
 
@@ -1038,7 +1057,11 @@ def recompute_component(self, inputs, block_variable, idx, prepared):
     if not self.forward_kwargs:
         dolfin_adjoint.backend.solve(lhs == rhs, func, bcs, solver_parameters={'linear_solver': 'mumps'})
     else:
-        dolfin_adjoint.backend.solve(lhs == rhs, func, bcs, **self.forward_kwargs)
+        if "krylov_solver_parameters" in self.forward_kwargs.keys():
+            solver_parameters={'linear_solver': self.forward_kwargs["krylov_method"], 'preconditioner': self.forward_kwargs["krylov_preconditioner"]}
+            dolfin_adjoint.backend.solve(lhs == rhs, func, bcs, solver_parameters=solver_parameters)
+        else:
+            dolfin_adjoint.backend.solve(lhs == rhs, func, bcs, **self.forward_kwargs)
         # print()
         # print("func = "+repr(np.mean(func.vector().get_local())))
 

windse/helper_functions.py

@@ -361,7 +361,10 @@ def dRForce(r,d_r): return (r*np.cos(np.pi*r)*(np.pi*d_r) + d_r*np.sin(np.pi*r))
 
 def UpdateActuatorLineForce(problem, mpi_u_fluid_constant, simTime_id, dt, turb_i, dfd=None, verbose=False):
 
+
+
     simTime = problem.simTime_list[simTime_id]
+    # print("debug data:", simTime,  mpi_u_fluid_constant.values())
 
     if verbose:
         print("Current Optimization Time: "+repr(simTime)+", Turbine #"+repr(turb_i))

windse_driver/driver.py

@@ -88,6 +88,8 @@ def run_action(params_loc=None):
     if params.rank == 0:
         print("Run Complete: {:1.2f} s".format(runtime))
 
+    params.comm.Barrier()
+
     return runtime
 
 def mesh_action(params_loc=None):