First prototype of the dune-precice module (dune-adapter) for preCICE

README.md

@@ -1,2 +1,34 @@
 # dune-adapter
 **experimental** preCICE-adapter for DUNE, a modular toolbox for solving partial differential equations
+
+## dune-precice
+A adapter module enabeling multiphysics simulation inside DUNE with the help of preCICE [1][2].
+
+The adapter can be build with:
+`<path-to-dune-common/bin/dunecontrol> --current all`
+
+For more detailed installation instructions have a look at:
+https://www.dune-project.org/doc/installation/
+
+## dune-precice-howto
+This module contains tutorial cases to get familiar with preCICE and the respective dune-adapter.
+
+The dune-adapter module depends on the DUNE core modules (version >= 2.7), dune-precice (precice version >= 2.0.0) and
+on dune-elastodynamics, dune-foamgrid (version >= 2.7) providing the necessary tools for structural simulation.
+
+The dune-elastodynamics module can be found here:
+https://github.com/maxfirmbach/dune-elastodynamics
+
+The executable for the perpendicular-flap tutorial case can be
+found in:
+`dune-precice/build-cmake/src/`
+after the module was build and needs to be copied to the solid-dune folder inside the respective
+tutorial case.
+
+<a id="1">[1]</a> 
+Firmbach M. (2021).
+[Aeroelastic simulation of slender wings for electric aircraft - A partitioned approach with DUNE and preCICE](https://mediatum.ub.tum.de/node?id=1609293), Master's Thesis
+
+<a id="2">[2]</a> 
+Firmbach M., Callies R. (2021).
+[Aeroelastic simulation of slender wings for electric aircraft - A partitioned approach with DUNE and preCICE](https://athene-forschung.unibw.de/138607), Conference Contribution

dune-precice-howto/CMakeLists.txt

@@ -0,0 +1,27 @@
+cmake_minimum_required(VERSION 3.13)
+project(dune-precice-howto CXX)
+
+if(NOT (dune-common_DIR OR dune-common_ROOT OR
+      "${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
+    string(REPLACE  ${PROJECT_NAME} dune-common dune-common_DIR
+      ${PROJECT_BINARY_DIR})
+endif()
+
+#find dune-common and set the module path
+find_package(dune-common REQUIRED)
+list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/modules"
+  ${dune-common_MODULE_PATH})
+
+#include the dune macros
+include(DuneMacros)
+
+# start a dune project with information from dune.module
+dune_project()
+
+dune_enable_all_packages()
+
+add_subdirectory(examples)
+add_subdirectory(cmake/modules)
+
+# finalize the dune project, e.g. generating config.h etc.
+finalize_dune_project(GENERATE_CONFIG_H_CMAKE)

dune-precice-howto/cmake/modules/CMakeLists.txt

@@ -0,0 +1,3 @@
+set(modules "DunePreciceHowtoMacros.cmake")
+
+install(FILES ${modules} DESTINATION ${DUNE_INSTALL_MODULEDIR})

dune-precice-howto/cmake/modules/DunePreciceHowtoMacros.cmake

@@ -0,0 +1 @@
+# File for module specific CMake tests.

dune-precice-howto/config.h.cmake

@@ -0,0 +1,45 @@
+/* begin dune-precice-howto
+   put the definitions for config.h specific to
+   your project here. Everything above will be
+   overwritten
+*/
+
+/* begin private */
+/* Name of package */
+#define PACKAGE "@DUNE_MOD_NAME@"
+
+/* Define to the address where bug reports for this package should be sent. */
+#define PACKAGE_BUGREPORT "@DUNE_MAINTAINER@"
+
+/* Define to the full name of this package. */
+#define PACKAGE_NAME "@DUNE_MOD_NAME@"
+
+/* Define to the full name and version of this package. */
+#define PACKAGE_STRING "@DUNE_MOD_NAME@ @DUNE_MOD_VERSION@"
+
+/* Define to the one symbol short name of this package. */
+#define PACKAGE_TARNAME "@DUNE_MOD_NAME@"
+
+/* Define to the home page for this package. */
+#define PACKAGE_URL "@DUNE_MOD_URL@"
+
+/* Define to the version of this package. */
+#define PACKAGE_VERSION "@DUNE_MOD_VERSION@"
+
+/* end private */
+
+/* Define to the version of dune-precice-howto */
+#define DUNE_PRECICE_HOWTO_VERSION "@DUNE_PRECICE_HOWTO_VERSION@"
+
+/* Define to the major version of dune-precice-howto */
+#define DUNE_PRECICE_HOWTO_VERSION_MAJOR @DUNE_PRECICE_HOWTO_VERSION_MAJOR@
+
+/* Define to the minor version of dune-precice-howto */
+#define DUNE_PRECICE_HOWTO_VERSION_MINOR @DUNE_PRECICE_HOWTO_VERSION_MINOR@
+
+/* Define to the revision of dune-precice-howto */
+#define DUNE_PRECICE_HOWTO_VERSION_REVISION @DUNE_PRECICE_HOWTO_VERSION_REVISION@
+
+/* end dune-precice-howto
+   Everything below here will be overwritten
+*/

dune-precice-howto/dune-precice-howto.pc.in

@@ -0,0 +1,15 @@
+prefix=@prefix@
+exec_prefix=@exec_prefix@
+libdir=@libdir@
+includedir=@includedir@
+CXX=@CXX@
+CC=@CC@
+DEPENDENCIES=@REQUIRES@
+
+Name: @PACKAGE_NAME@
+Version: @VERSION@
+Description: dune-precice-howto module
+URL: http://dune-project.org/
+Requires: dune-common dune-geometry dune-uggrid dune-grid dune-istl dune-functions dune-elastodynamics dune-precice
+Libs: -L${libdir}
+Cflags: -I${includedir}

dune-precice-howto/dune.module

@@ -0,0 +1,12 @@
+################################
+# Dune module information file #
+################################
+
+# Name of the module
+Module: dune-precice-howto
+Version: 1.0
+Maintainer: [email protected]
+# Required build dependencies
+Depends: dune-common dune-geometry dune-uggrid dune-grid dune-istl dune-functions dune-elastodynamics dune-precice
+# Optional build dependencies
+#Suggests:

dune-precice-howto/examples/CMakeLists.txt

@@ -0,0 +1,12 @@
+find_package(precice REQUIRED CONFIG)
+
+set(programs
+		dune-perpendicular-flap)
+
+#include precice
+find_package(precice)
+
+foreach(_program ${programs})
+	add_executable(${_program} ${_program}.cc)
+	target_link_libraries(${_program} PRIVATE precice::precice)
+endforeach()

dune-precice-howto/examples/dune-perpendicular-flap.cc

@@ -0,0 +1,222 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#include <config.h>
+
+#include <stdio.h>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include <dune/istl/matrixmarket.hh>
+#include <dune/istl/matrix.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bdmatrix.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/owneroverlapcopy.hh>
+
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/boundarydofs.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
+#include <dune/elastodynamics/assemblers/operatorassembler.hh>
+#include <dune/elastodynamics/assemblers/stiffnessassembler.hh>
+#include <dune/elastodynamics/assemblers/consistentmassassembler.hh>
+
+#include <dune/elastodynamics/timesteppers/coefficients.hh>
+#include <dune/elastodynamics/timesteppers/timestepcontroller.hh>
+#include <dune/elastodynamics/timesteppers/newmark.hh>
+
+#include <dune-precice/couplinginterface.hh>
+
+using namespace Dune;
+
+struct couplingParameters {
+
+  const std::string config_file      = "../precice-config.xml";
+  const std::string participant_name = "Solid";
+  const std::string mesh_name        = "Solid-Mesh";
+  const std::string read_data_name   = "Force";
+  const std::string write_data_name  = "Displacement";
+
+};
+
+int main(int argc, char** argv) {
+
+  // set up MPI
+  const MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
+  int mpiRank = mpiHelper.rank();
+  int mpiSize = mpiHelper.size();
+
+  const int dim = 2;
+  const int p = 2;
+
+  // load and generate grid
+  using Grid = UGGrid<dim>;
+  using GridView = Grid::LeafGridView;
+
+  GridFactory<Grid> factory;
+  GmshReader<Grid>::read(factory, "Solid.msh");
+  std::shared_ptr<Grid> grid(factory.createGrid());    
+
+  // partition grid on different processors
+  GridView gridView = grid->leafGridView();
+  std::cout << "\nThere are " << grid->leafGridView().comm().size() << " processes active!" << std::endl;
+
+  // generate Basis
+  using namespace Functions::BasisFactory;
+  auto basis = makeBasis(gridView, power<dim>(lagrange<p>()));
+  using Basis = decltype(basis);
+
+  // define operators needed
+  using operatorType = BCRSMatrix<FieldMatrix<double, dim, dim>>;
+  using blockVector  = BlockVector<FieldVector<double, dim>>;
+
+  // assemble problem
+  Elastodynamics::OperatorAssembler<Basis> operatorAssembler(basis);
+
+  double E = 4000000.0, nu = 0.3;
+  operatorType stiffnessMatrix;
+  operatorAssembler.initialize(stiffnessMatrix);
+  Elastodynamics::StiffnessAssembler stiffnessAssembler(E, nu);
+  operatorAssembler.assemble(stiffnessAssembler, stiffnessMatrix, false);
+
+  double rho = 3000.0;
+  operatorType massMatrix;
+  operatorAssembler.initialize(massMatrix);
+  Elastodynamics::ConsistentMassAssembler massAssembler(rho);
+  operatorAssembler.assemble(massAssembler, massMatrix, false);
+
+  blockVector loadVector(basis.size()), displacementVector(basis.size()), velocityVector(basis.size()), accelerationVector(basis.size());
+  loadVector = 0.0, displacementVector = 0.0, velocityVector = 0.0, accelerationVector = 0.0;
+
+  // state quantatiy vector
+  std::vector<blockVector*> stateQuantaties = {&displacementVector, &velocityVector, &accelerationVector};
+
+  // set dirichlet boundary
+  BitSetVector<dim> dirichletDofs(basis.size(), false);
+  auto dirichletPredicate = [](auto position) {return position[1] < 0.00001;};  
+  Functions::interpolate(basis, dirichletDofs, dirichletPredicate);
+
+  DiagonalMatrix<double, dim> I(1.0);
+  FieldMatrix<double, dim> OM(0.0);
+
+  // set fixed Dirichlet entries in stiffnessmatrix
+  for( int i=0; i<stiffnessMatrix.N(); i++) {
+    if(dirichletDofs[i][0]) {
+      auto cIt = stiffnessMatrix[i].begin();
+      auto cEndIt = stiffnessMatrix[i].end();  
+      for( ; cIt!=cEndIt; ++cIt) {
+        *cIt = (i==cIt.index()) ? I : OM;
+      }
+    }
+  }
+
+  // set fixed Dirichlet entries in massmatrix
+  for( int i=0; i<massMatrix.N(); i++) {
+    if(dirichletDofs[i][0]) {
+      auto cIt = massMatrix[i].begin();
+      auto cEndIt = massMatrix[i].end();  
+      for( ; cIt!=cEndIt; ++cIt) {
+        *cIt = (i==cIt.index()) ? I : OM;
+      }
+    }
+  }
+
+  FieldVector<double, dim> OV(0.0);
+
+  // set fixed Dirichlet entries in load vector
+  for( int i=0; i<loadVector.N(); i++) {
+    if(dirichletDofs[i][0]) {
+      loadVector[i] = OV;
+    }
+  }
+
+  // set neumann boundary
+  BlockVector<FieldVector<bool, dim>> neumannDofs(basis.size());
+  neumannDofs = false;
+  auto neumannPredicate = [](auto position) {return position[1] > -3.000001 || position[0] < 2.955555 || position[0] > 3.055555;};  
+
+  Functions::interpolate(basis, neumannDofs, neumannPredicate);
+
+  // generate global displacement function
+  using displacementRange = Dune::FieldVector<double, dim>;
+  auto displacementFunction = Functions::makeDiscreteGlobalBasisFunction<displacementRange> (basis, displacementVector);
+
+  // generate output writer
+  double dt_out = 0.1;
+  auto vtkWriter = std::make_shared<SubsamplingVTKWriter<GridView>> (gridView, refinementLevels(0));
+  VTKSequenceWriter<GridView> vtkSequenceWriter(vtkWriter, "solid");
+  vtkWriter->addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
+  vtkSequenceWriter.write(0.0);
+
+  // setup looping
+  double t = 0.0;
+  double dt = 0.01;
+  double precice_dt = 0.0;
+  int iteration_count = 0;
+
+  // set up Newmark method
+  FixedStepController fixed(t, dt);
+  NewmarkCoefficients coefficients = ConstantAcceleration();
+  Newmark<operatorType, blockVector> newmark(massMatrix, stiffnessMatrix, coefficients, fixed);
+  newmark.initialize(accelerationVector, loadVector);
+
+  // set up coupling interface
+  couplingParameters parameters; 
+  preCICE::CouplingInterface<dim, blockVector, couplingParameters> couplingInterface(parameters, neumannDofs, mpiRank, mpiSize); 
+  precice_dt = couplingInterface.initialize(basis);
+
+  while(couplingInterface.is_coupling_ongoing()) {
+
+    if(couplingInterface.is_save_required()) {
+        couplingInterface.save_current_state(stateQuantaties, t, iteration_count);
+        couplingInterface.mark_save_fullfilled();
+    }
+
+    // get neumann values from fluid solver
+    couplingInterface.read_blockvector_data(loadVector);
+
+    // setup fixed dirichlet values
+    for( int i=0; i<loadVector.N(); i++) {
+      if(dirichletDofs[i][0]) {
+        loadVector[i] = OV;
+      }
+    }
+
+    newmark.step(displacementVector, velocityVector, accelerationVector, loadVector);
+
+    dt = std::min(precice_dt, dt);
+
+    couplingInterface.write_blockvector_data(displacementVector);
+    precice_dt = couplingInterface.advance(dt);
+
+    if(couplingInterface.is_load_required()) {
+        couplingInterface.reload_old_state(stateQuantaties, t, iteration_count);
+        couplingInterface.mark_load_fullfilled();
+    }
+    else {
+      t += dt;
+      ++iteration_count;
+      std::cout << "=== iteration: " << iteration_count << " time: " << t << std::endl;
+    }
+
+    if(couplingInterface.is_time_window_complete()) {
+      double tol = 10e-5;
+      if(std::abs(std::fmod((t+tol), dt_out)) < 2.0*tol) {
+        std::cout << "output vtk for time = " << t << "\n" << std::endl;
+        vtkSequenceWriter.write(t);
+      }
+    }
+  }
+
+  couplingInterface.finalize();
+
+}