Merge pull request #18 from KVSlab/update-tests

Update tests [WIP]

Author: hkjeldsberg

.github/workflows/check_and_test_package.yml

@@ -87,7 +87,7 @@ jobs:
         if: steps.cache.outputs.cache-hit != 'true'
 
       - name: Install OasisMove
-        run: python3 -m pip install .[test]
+        run: python3 -m pip install --editable .[test]
 
       - name: Run tests
         run: python3 -m pytest tests

src/oasismove/problems/NSfracStep/DrivenCavity.py

@@ -6,7 +6,7 @@
 
 
 # Override some problem specific parameters
-def problem_parameters(NS_parameters, **NS_namespace):
+def problem_parameters(NS_parameters, scalar_components, **NS_namespace):
     NS_parameters.update(
         # Mesh parameters
         Nx=50,
@@ -28,6 +28,10 @@ def problem_parameters(NS_parameters, **NS_namespace):
         pressure_degree=1,
         use_krylov_solvers=True)
 
+    scalar_components += ["alfa", "beta"]
+    Schmidt["alfa"] = 1.
+    Schmidt["beta"] = 10.
+
 
 # Create a mesh
 def mesh(Nx=50, Ny=50, **params):
@@ -39,12 +43,16 @@ def mesh(Nx=50, Ny=50, **params):
 def create_bcs(V, **NS_namespace):
     noslip = "std::abs(x[0]*x[1]*(1-x[0]))<1e-8"
     top = "std::abs(x[1]-1) < 1e-8"
+    bottom = "std::abs(x[1]) < 1e-8"
     bc0 = DirichletBC(V, 0, noslip)
     bc00 = DirichletBC(V, 1, top)
     bc01 = DirichletBC(V, 0, top)
+    bcbeta = DirichletBC(V, 1, bottom)
     return dict(u0=[bc00, bc0],
                 u1=[bc01, bc0],
-                p=[])
+                p=[],
+                alfa=[bc00],
+                beta=[bcbeta])
 
 
 def initialize(x_1, x_2, bcs, **NS_namespace):

src/oasismove/problems/NSfracStep/MovingCylinder.py

@@ -166,8 +166,8 @@ def update_boundary_conditions(t, NS_expressions, **NS_namespace):
     NS_expressions["circle_y"].t = t
 
 
-def temporal_hook(t, St, F, A_ratio, tstep, save_solution_frequency, forces, q_, u_inf, D, viz_u, viz_p, u_vec,
-                  newfolder, p_, u_, **NS_namespace):
+def temporal_hook(mesh, t, dt, nu, St, F, A_ratio, tstep, save_solution_frequency, forces, q_, u_inf, D, viz_u, viz_p,
+                  u_vec, newfolder, p_, u_, **NS_namespace):
     # Save fluid velocity and pressure solution
     if tstep % save_solution_frequency == 0:
         assign(u_vec.sub(0), u_[0])
@@ -208,3 +208,11 @@ def temporal_hook(t, St, F, A_ratio, tstep, save_solution_frequency, forces, q_,
         lift_coeff = sum(drag_and_lift[1::2]) / factor
         data = [t, tstep, drag_coeff, lift_coeff, pressure_coeff]
         write_data_to_file(newfolder, data, "forces.txt")
+
+    # Compute mean velocity and Reynolds number at inlet
+    if tstep % 10 == 0:
+        h = mesh.hmin()
+        L = 62 * D
+        compute_flow_quantities(u_, L, nu, mesh, t, tstep, dt, h, outlet_area=1, boundary=None, outlet_ids=[],
+                                inlet_ids=[], id_wall=0, period=1.0, newfolder=None, dynamic_mesh=False,
+                                write_to_file=False)

src/oasismove/problems/NSfracStep/MovingTaylorGreen3D.py

@@ -4,6 +4,8 @@
 from oasismove.problems.NSfracStep import *
 from oasismove.problems.NSfracStep.MovingCommon import get_visualization_writers
 
+comm = MPI.comm_world
+
 
 def problem_parameters(commandline_kwargs, NS_parameters, NS_expressions, **NS_namespace):
     """
@@ -175,10 +177,17 @@ def update_boundary_conditions(t, NS_expressions, **NS_namespace):
             value.t = t
 
 
-def temporal_hook(u_, save_solution_frequency, u_vec, viz_u, tstep, t, **NS_namespace):
+def temporal_hook(nu, L, dt, mesh, u_, save_solution_frequency, u_vec, viz_u, tstep, t, **NS_namespace):
     # Save velocity and pressure
     if tstep % save_solution_frequency == 0:
         for i in range(3):
             assign(u_vec.sub(i), u_[i])
 
         viz_u.write(u_vec, t)
+
+    # Compute mean velocity and Reynolds number at inlet
+    if tstep % 10 == 0:
+        h = mesh.hmin()
+        compute_flow_quantities(u_, L, nu, mesh, t, tstep, dt, h, outlet_area=1, boundary=None, outlet_ids=[],
+                                inlet_ids=[], id_wall=0, period=1.0, newfolder=None, dynamic_mesh=False,
+                                write_to_file=False)

src/oasismove/problems/NSfracStep/MovingWall.py

@@ -1,6 +1,8 @@
 from oasismove.problems.NSfracStep import *
 from oasismove.problems.NSfracStep.MovingCommon import get_coordinate_map, get_visualization_writers
 
+comm = MPI.comm_world
+
 
 def problem_parameters(NS_parameters, **NS_namespace):
     """
@@ -196,11 +198,20 @@ def update_boundary_conditions(t, NS_expressions, **NS_namespace):
             value.update(t)
 
 
-def temporal_hook(t, tstep, save_solution_frequency, p_, u_, viz_u, viz_p, u_vec, **NS_namespace):
+def temporal_hook(mesh, dt, h0, eps, nu, t, tstep, save_solution_frequency, p_, u_, viz_u, viz_p, u_vec,
+                  **NS_namespace):
     # Write solution at time t
     if tstep % save_solution_frequency == 0:
         assign(u_vec.sub(0), u_[0])
         assign(u_vec.sub(1), u_[1])
 
         viz_u.write(u_vec, t)
         viz_p.write(p_, t)
+
+    # Compute mean velocity and Reynolds number at inlet
+    if tstep % 10 == 0:
+        h = mesh.hmin()
+        L = h0 * (1 + eps)
+        compute_flow_quantities(u_, L, nu, mesh, t, tstep, dt, h, outlet_area=1, boundary=None, outlet_ids=[],
+                                inlet_ids=[], id_wall=0, period=1.0, newfolder=None, dynamic_mesh=False,
+                                write_to_file=False)

src/oasismove/problems/NSfracStep/probe_v2/Probe.cpp

@@ -0,0 +1,186 @@
+#include "Probe.h"
+
+using namespace dolfin;
+
+Probe::~Probe() 
+{
+  clear(); 
+}
+
+Probe::Probe(const Array<double>& x, const FunctionSpace& V) :
+  _element(V.element()), _num_evals(0)
+{
+  auto mesh = V.mesh();
+  std::size_t gdim = mesh->geometry().dim();
+  _x.resize(3);
+
+  // Find the cell that contains probe
+  const Point point(gdim, x.data());
+  unsigned int cell_id = mesh->bounding_box_tree()->compute_first_entity_collision(point);
+
+  // If the cell is on this process, then create an instance 
+  // of the Probe class. Otherwise raise a dolfin_error.
+  if (cell_id != std::numeric_limits<unsigned int>::max())
+  {
+    // Store position of probe
+    for (std::size_t i = 0; i < 3; i++) 
+      _x[i] = (i < gdim ? x[i] : 0.0);
+
+    // Compute in tensor (one for scalar function, . . .)
+    _value_size_loc = 1;
+    for (std::size_t i = 0; i < _element->value_rank(); i++)
+       _value_size_loc *= _element->value_dimension(i);
+
+    _probes.resize(_value_size_loc);
+
+    // Create cell that contains point
+    dolfin_cell.reset(new Cell(*mesh, cell_id));
+    dolfin_cell->get_cell_data(ufc_cell);
+    
+    coefficients.resize(_element->space_dimension());
+    
+    // Cell vertices
+    dolfin_cell->get_vertex_coordinates(vertex_coordinates);
+        
+    // Create work vector for basis
+    basis_matrix.resize(_value_size_loc);
+    for (std::size_t i = 0; i < _value_size_loc; ++i)
+      basis_matrix[i].resize(_element->space_dimension());
+        
+    std::vector<double> basis(_value_size_loc);
+    const int cell_orientation = 0;
+    for (std::size_t i = 0; i < _element->space_dimension(); ++i)
+    {
+      _element->evaluate_basis(i, &basis[0], &x[0], 
+                               vertex_coordinates.data(), 
+                               cell_orientation);
+      for (std::size_t j = 0; j < _value_size_loc; ++j)
+        basis_matrix[j][i] = basis[j];
+    }
+  }  
+  else
+  {  
+    dolfin_error("Probe.cpp","set probe","Probe is not found on processor");
+  }
+}
+//
+Probe::Probe(const Probe& p)
+{
+  basis_matrix = p.basis_matrix;
+  coefficients = p.coefficients;
+  _x = p._x;
+  vertex_coordinates = p.vertex_coordinates;
+  _element = p._element;
+  dolfin_cell.reset(new Cell(p.dolfin_cell->mesh(), p.dolfin_cell->index()));
+  ufc_cell = p.ufc_cell;
+  _value_size_loc = p._value_size_loc;
+  _num_evals = p._num_evals;
+  _probes = p._probes;
+}
+//
+void Probe::eval(const Function& u)
+{
+  // Restrict function to cell
+  u.restrict(&coefficients[0], *_element, *dolfin_cell,
+              vertex_coordinates.data(), ufc_cell);
+
+  // Make room for one more evaluation
+  for (std::size_t j = 0; j < _value_size_loc; j++)
+    _probes[j].push_back(0.);
+  
+  // Compute linear combination
+  for (std::size_t i = 0; i < _element->space_dimension(); i++)
+  {
+    for (std::size_t j = 0; j < _value_size_loc; j++)
+      _probes[j][_num_evals] += coefficients[i]*basis_matrix[j][i];
+  }
+  
+  _num_evals++;
+}
+// Remove one instance of the probe
+void Probe::erase_snapshot(std::size_t i)
+{
+  for (std::size_t j = 0; j < _value_size_loc; j++)
+    _probes[j].erase(_probes[j].begin()+i);  
+  
+  _num_evals--;
+}
+// Reset probe by removing all values
+void Probe::clear()
+{
+  for (std::size_t j = 0; j < _value_size_loc; j++)
+    _probes[j].clear();  
+  
+  _num_evals = 0;
+}
+// Return probe values for chosen component of tensor
+std::vector<double> Probe::get_probe_sub(std::size_t i)
+{
+  return _probes[i];
+}
+// Return probe values for entire tensor at snapshot
+std::vector<double> Probe::get_probe_at_snapshot(std::size_t i)
+{
+  std::vector<double> x(_value_size_loc);
+  for (std::size_t j = 0; j < _value_size_loc; j++)
+    x[j] = _probes[j][i];
+  return x;
+}
+// Return robe value for given component at given snapshot
+double Probe::get_probe_component_and_snapshot(std::size_t comp, std::size_t i)
+{
+  return _probes[comp][i];
+}
+// Return coordinates of probe
+std::vector<double> Probe::coordinates()
+{
+  std::vector<double> x(_x);
+  return x;
+}
+//
+void Probe::dump(std::size_t i, std::string filename, std::size_t id)
+{
+  //std::ofstream fp(filename.c_str(), std::ios_base::app);
+  std::ofstream fp(filename.c_str(), std::ios_base::out);
+  std::ostringstream ss;
+  ss << std::scientific;
+  ss << "Probe id = " << id << std::endl;
+  ss << "Number of evaluations = " << number_of_evaluations() << std::endl;
+  ss << "Coordinates:" << std::endl;
+  ss << _x[0] << " " << _x[1] << " " << _x[2] << std::endl; 
+  ss << "Values for component " << i << std::endl;
+  for (std::size_t j = 0; j < _probes[i].size(); j++)
+      ss << _probes[i][j] << std::endl;
+  ss << std::endl;
+  cout << ss.str() << endl;
+  fp << ss.str();
+}
+void Probe::dump(std::string filename, std::size_t id)
+{
+  //std::ofstream fp(filename.c_str(), std::ios_base::app);
+  std::ofstream fp(filename.c_str(), std::ios_base::out);
+  std::ostringstream ss;
+  ss << std::scientific;
+  ss << "Probe id = " << id << std::endl;
+  ss << "Number of evaluations = " << number_of_evaluations() << std::endl;
+  ss << "Coordinates:" << std::endl;
+  ss << _x[0] << " " << _x[1] << " " << _x[2] << std::endl; 
+  ss << "Values for all components:" << std::endl;    
+  std::size_t N = _probes[0].size();
+  for (std::size_t j = 0; j < N; j++)
+  {
+    for (std::size_t i = 0; i < _value_size_loc; i++)
+      ss << _probes[i][j] << " " ;
+    ss << std::endl;
+  }
+  cout << ss.str() << endl;
+  fp << ss.str();
+}
+//
+void Probe::restart_probe(const Array<double>& u)
+{
+  for (std::size_t j = 0; j < _value_size_loc; j++)
+    _probes[j].push_back(u[j]);
+}
+
+