Implement svFSI using C++

[ktbolt]

While Fortran is still used for scientific and engineering problems we feel that it would be better to implement the FSI solver using C++ to provide a more extensible framework (e.g. plugins for custom BCs) and better capabilities for text file processing (e.g. XML).

Although Fortran is easier to learn than C++ and has nice array support most students don't have any experience programming in Fortran.

We will perform the conversion in a series of steps. The first step is to implement a C++ framework providing a Fortran-callable interface used to store and retrieve data. Each Fortran module or other independent code sections will then be incrementally converted to C++.

incrementally converting sections of independent code.

[ktbolt]

I've written some Fortran code to prototype ways to interface the svFSI Fortran code to a c++ framework.

I've implemented a Fortran simulation_interface module as an interface to the c++ Simulation class used to read an FSI XML file and VTK mesh files.

The code below shows how to use the simulation_interface module to create a Simulation object, use it to read in an svFSI XML input parameters file, load the mesh (VTK files) and get an array of nodal coordinates for that mesh

! Test interfacing to the Simulation c++ class.
! 
program test_sim_interface

  use simulation_interface
  implicit none
  type(simulation) :: sim_interface       
  real(c_double), pointer :: node_coords(:,:)

  ! Get the solver input xml file name.
  character (len=32) :: in_file_name
  call getarg(1, in_file_name)

  ! Create a c++ Simulation object.
  sim_interface = simulation()

  ! Read in the solver input xml file.
  call sim_interface % read_parameters(in_file_name)

  ! Load the mesh.
  call sim_interface % load_mesh()

  ! Get node coordinates. 
  call sim_interface % get_node_coords(node_coords)

  num_coords = ubound(node_coords,1)
  do i = 1, num_coords
      print *, "[main] node_coords(1:): ", node_coords(i,1), node_coords(i,2), node_coords(i,3)
  end do

The node_coords array can then be used in the svFSI code as it normally would if it was defined in the current svFSI code using

!   Position coordinates
REAL(KIND=RKIND), ALLOCATABLE :: x(:,:)

[ktbolt]

I've integrated the C++ prototype code into svFSI, creating a C++ Simulation object and read in a solver XML file

      PROGRAM MAIN

      USE COMMOD
      USE ALLFUN
      use compare_cpp_fortran

      IMPLICIT NONE
.
.
.
      ! Read in the solver input XML file.
      !
      ! Assume that the XML file has the same name as the .txt 
      ! input solver file.
      !
      call getarg(1, in_file_name)
      str_len = len(trim(in_file_name))
      write (*,*) '[main] Len input file name: ', str_len
      in_file_name = trim(in_file_name(1:str_len-3)) // 'xml'
      write (*,*) '[main] Input XML file name: ', in_file_name
      call sim_interface % read_parameters(in_file_name)

      ! Load the mesh.
      call sim_interface % load_mesh()

      ! Create auxillary data needed for a simulation. 
      call sim_interface % create_aux_data()

      ! Read in solver input file, storing parameter values in 'COMMOD'. 
      !
      print *, "[MAIN] Call READFILES -->"
 101  CALL READFILES
      print *, "[MAIN] --> READFILES return"
      print *, "[MAIN] nsd: ", nsd

      call compare_mesh()

I've also added a compare_cpp_fortran module used to compare Fortran and C++ data.

[vvedula22]

@ktbolt I merged Chi's P2P1 stuff for fluid and FSI. There were some other modifications in which I extended P2P1 for non-Newtonian flows and weakly applied Dirichlet BCs. Please use the commit https://github.com/SimVascular/svFSI/commit/e258b9fed7f3f2745bdc5d089c8a37fe1e1fd877 for your C++ conversion. I have tested this commit on all fluid and FSI cases in the svFSI-Tests repository. All the code additions we do subsequently including Martin's G&R stuff, any additions to cardiac mechanics, shells, IB, etc. could be converted later.

[ktbolt]

@vvedula22 Great! I will get back to the conversion.

[ktbolt]

I've redone the XML file format to mimic the layout and names currently used by svFSI

<GeneralSimulationParameters>
  <Continue_previous_simulation> 0 </Continue_previous_simulation>
  <Number_of_spatial_dimensions> 3 </Number_of_spatial_dimensions>
  <Number_of_time_steps> 800 </Number_of_time_steps>
  <Warning> 0 </Warning>
  <Debug> 0 </Debug>
</GeneralSimulationParameters>

<Add_mesh name="msh" >
  <Mesh_file_path> mesh-complete/mesh-complete.mesh.vtu </Mesh_file_path>
  <Add_face name="lumen_outlet">
      <Face_file_path> mesh-complete/mesh-surfaces/lumen_outlet.vtp </Face_file_path>
  </Add_face>
  <Add_face name="lumen_wall">
      <Face_file_path> mesh-complete/mesh-surfaces/lumen_wall.vtp </Face_file_path>
  </Add_face>
</Add_mesh>

<Add_equation name="fluid" >
   <Density> 1.06 </Density>
   <Viscosity model="Constant" >
     <Value> 0.04 </Value>
   </Viscosity>
   <Coupled> true </Coupled>
   <Min_iterations> 3 </Min_iterations>
.
.
.
</Add_equation>

I also rewrote the XML reader to remove the jumble of parameter names and such.

[ktbolt]

I've created an Implement-svFSI-using-cpp_19 branch from my svFSI fork, integrated new the XML reader and added the previous C++ interface code.

I'm now creating C++ classes that duplicate Fortran modules, retain the same names and such to make it easy to match the Fortran data. I also created a simple C++ Array class to store matrices.

I've converted the types in theCOMMOD module, the mother of all modules, into C++ classes. For example

TYPE fsType
       LOGICAL lShpF
       INTEGER(KIND=IKIND) :: eType = eType_NA
       INTEGER(KIND=IKIND) eNoN
...
       REAL(KIND=RKIND), ALLOCATABLE :: N(:,:)
       REAL(KIND=RKIND), ALLOCATABLE :: Nb(:,:)
END TYPE fsType

is converted to

class fsType {
  bool lShpF;
  int eType;
  int eNoN;
...
  Array<double> N;
  Array<double> Nb;
};

I've stated to duplicate the Fortran subroutines to maintain a similar flow of control and calling sequence.

[ktbolt]

I've added new XML elements to the parser to support fsi domains and more complicated BCs

<Add_equation name="FSI" >
   <Coupled> true </Coupled>
   <Min_iterations> 3 </Min_iterations>
   <Max_iterations> 10 </Max_iterations>
   <Tolerance> 1e-3 </Tolerance>

   <Domain id="0" >
       <Equation> fluid </Equation>
       <Density> 1.0 </Density>
       <Viscosity model="Constant" >
           <Value> 0.04 </Value>
       </Viscosity>
       <Backflow_stabilization_coefficient> 0.2 </Backflow_stabilization_coefficient>
   </Domain>

   <Add_BC name="wall_inlet" >
      <Type> Dir </Type>
      <Value> 0.0 </Value>
      <Impose_on_state_variable_integral> true </Impose_on_state_variable_integral>
      <Zero_out_perimeter> false </Zero_out_perimeter>
      <Effective_direction> (0, 0, 1) </Effective_direction>
   </Add_BC>

I've added returning mesh data from the C++ to Fortran for comparison, nodal coordinates match fine but element connectivity was wrong, elements ordered differently, missed important step later where Reordering element connectivity occurred and node IDs were incremented by 1 but not in the same place.

So now I've reorganized the C++ code to better replicate the Fortran organization and flow of control, replicating where Fortran subroutines are called from and their name.

[ktbolt]

I've reorganized the code and have replicated the Fortran code I've come across so far as I track the control flow for a CFD simulation.

consts.h  consts.cpp
load_msh.h  load_msh.cpp
nn.h. nn.cpp. nn_elem_gip.h. nn_elem_gnn.h  nn_elem_nn_bnds.h. nn_elem_props.h
read_files.h  read_files.cpp
read_msh.h. read_msh.cpp
vtk_xml.h  vtk_xml.cpp
vtk_xml_parser.h  vtk_xml_parser.cpp

It's very easy to get mixed up here so I'm documenting all of this in https://github.com/ktbolt/svFSI/blob/Implement-svFSI-using-cpp_19/Code/Source/svFSI_cinterface/README.md.

Even though a lot of the element code is screaming object orient me I'm mostly following the Fortran implementation. I did remove a lot of element type and element nodes case statements, such as SELECT CASE (lM%eNoN) and SELECT CASE (lFa%eNoN), that clutters up the the code, replaced them with maps of lambda functions. All of this type of thing will be replaced with element classes later on.

[ktbolt]

I've finished converting all the Fortran code called from LOADMSH: reading meshes and faces, checking and reordering element connectivity, setting Gauss integration and shape function element data.

Subtracting 1 from array and vector indices, etc. is totally error prone of course but have Array and Vector class checks to track mistakes down.

I've been adding Array and Vector class operators as needed to try to replicate some of what Fortran supports. For example, the Fortran code

xl = lM % x(:, lM % gIEN(:,e))
v(:,1) = xl(:,2) - xl(:,1)
v(:,2) = xl(:,3) - xl(:,2)
v(:,3) = xl(:,4) - xl(:,3)
v(:,4) = CROSS(v(:,1:2))
sn(1)  = SGN(SUM(v(:,3)*v(:,4)))

is similarly implemented like this in C++

auto elem_conn = mesh.gIEN.col(e);     
xl = mesh.x.get_cols(elem_conn);
v1 = xl[1] - xl[0];
v2 = xl[2] - xl[1];
v3 = xl[3] - xl[2];
auto v4 = v1.cross(v2);
int sn = utils::sign(v3.dot(v4));

And the award for the most obscure Fortran subroutine name (so far) goes to READENDNLFF.

[vvedula22]

And the award for the most obscure Fortran subroutine name (so far) goes to READENDNLFF.

Ha ha! this refers to reading end (or boundary) nodes of a 1D fiber mesh from a file generated using the SV Purkinje plugin :) You will see many more later :)

[ktbolt]

I've finished converting the SETPROJECTOR(), FINDFACE(), MATCHFACES(), PULLSTACK() and PUSHSTACK() functions. This code is tricky because of the use of integer indexes into arrays storing node IDs , array bounds from [0,N-1] and such.

Need to be careful with lines like this in Fortran

stk%n = stk%n + 1
stk%v(stk%n) = iVal

to get the indexing right in C++.

I've added definitions to the Vector class to support STL iterator operations, for example for a Vector values you can do

// Iterate over values.
for (auto v : values) {
}

// Search for a value.
std::find(values.begin(), values.end(), value) == values.end());

I've replaced using std::vector in Array and other code with Vector to provide more uniform operations. For example

Vector<Vector<T>> get_cols(const Vector<int>& columns) const

Vector<T> rows(const int row, const Vector<int>& indexes) const

Note that some Fortran functions like PULLSTACK() and PUSHSTACK() can be replaced with the std::stack container. I've implemented them using Vector for now to replicate what's being done in Fortran.

[ktbolt]

I've finished converting the 500 line SUBROUTINE READMSH and the 12 subroutines that it calls to read in all sorts of things, added a bunch of parameters for prestress, initial field values, fiber directions, etc.

I've converted solver .inp to .xml files and tested on the following svFSI-Tests

07-fsi/ale/03-pipe_3D/
07-fsi/ale/04-pipe_prestress/02-deform-ale/02.2-fsi-ALE
06-ustruct/03-LV-Guccione-active
04-fluid/01-pipe3D_RCR/

Fortran and C++ codes data seem to match so far.

[ktbolt]

I'm converting READEQ now, also converting svFSI_GenBC.inp to XML for testing, will need to add new parameter for Couple_to_genBC.

And it finally dawned on me what the l means in variable names like lEq, lM, etc. ... local! Am I slow?

[ktbolt]

I've converted READEQ, READDOMAIN, READLS, READCEP, READMATMODEL, READVISCMODEL and FFT.

I've also added a bunch of enum classes to replicate the Fortran PARAMETER definitions in CONSTS.f. The Fortran

      INTEGER(KIND=IKIND), PARAMETER :: prop_NA = 0, fluid_density = 1,
     2   solid_density = 2, solid_viscosity = 3, elasticity_modulus = 4,
     3   poisson_ratio = 5, conductivity = 6, f_x = 7, f_y = 8, f_z = 9,
     4   backflow_stab = 10, source_term = 11, damping = 12,
     5   shell_thickness = 13, ctau_M = 14, ctau_C = 15

is implement in C++ like this

enum class PhysicalProperyTypes
{
  NA = 0,
  fluid_density = 1,
  solid_density = 2,
  solid_viscosity = 3,
  elasticity_modulus = 4,
  poisson_ratio = 5,
  conductivity = 6,
  f_x = 7,
  f_y = 8,
  f_z = 9,
  backflow_stab = 10,
  source_term = 11,
  damping = 12,
  shell_thickness = 13,
  ctau_M = 14,
  ctau_C = 15
};

These enums are referenced using PhysicalProperyTypes::fluid_density for example.

I've also converted more solver .inp files to .xml for testing.

[ktbolt]

@vvedula22 I'm testing setting CMM parameters using https://github.com/SimVascular/svFSI-Tests/blob/master/07-fsi/cmm/01-pipe_RCR/02-deform-wall_inflate/01-inflate/svFSI.inp.

The mesh is defined using

Add mesh: wall {
   Set mesh as shell: t
   Mesh file path: ../../mesh/walls_combined.vtp
}

The Fortran code is calling SUBROUTINE READVTU(lM, fName) to read in the VTP file

     CALL loadVTK(vtu, fName, iStat)
     IF (iStat .LT. 0) err = "VTU file read error (init)"

This is confusing because the name and implementation of READVTU implies that it should only read and process VTU files but it seems to be reading and processing VTP files as well. What's going on here?

[vvedula22]

@ktbolt This is a special case because we are loading the mesh as a shell surface where the triangulated wall surface is treated as a svFSI mesh and not a face. Likewise, the svFSI face for the shell surface would be the edges of the shell. Ideally, while creating the example, I should have copied the walls_combined file to a .vtu extension to maintain consistency. However, I didn't want to recreate files when not needed.

Regarding READVTU's ability to process VTP files, it is certainly possible because the Fortran VTK library is used to process only basic information in the vtu/vtp files such as position coordinates, connectivity, and any field variables. We don't do any other VTK operation that is specific to the data type, e.g., computing normals, extracting faces, etc.

Internally, the vtkXMLMod still loads the walls_combined file as .vtp file and looks for all the VTP commands such as NumberOfPolys instead of NumberOfCells, etc. However, READVTU only extracts coordinates and connectivity from this structure. It will not access other VTP data such as GlobalElementID.

To summarize, one can load a VTP file using READVTU as long as you don't need the additional vtp-specific field variables such as GlobalElementID or GlobalNodeID. Likewise, one can also load a VTU file using READVTP provided you have these additional field variables in the .vtu file and are required and appropriately processed for the problem.

[ktbolt]

@vvedula22 Thanks for the detailed reply!

I implemented READVTU to just read .vtu files, will think about redoing that or better yet I will convert the walls_combined.vtp file to walls_combined.vtu for my testing.

It's not clear why SV even uses vtp files to store mesh data, just overcomplicates processing mesh files, should use vtu for all mesh data.

[ktbolt]

I've converted the 500 line PARTMSH Fortran function that partitions and distributes a mesh across a set of processors, a most brutal experience, a 0-based Fortran array, figuring out which variables are counters and which are element IDs, redoing the MPI calls, and converting the CMMOD module implemented in COMU.f.

It will be interesting to reimplement this code not just following the Fortran implementation but using C++ containers, maybe can simplify this code. Just using descriptive variable names and code decomposition would help a lot.

[ktbolt]

@vvedula22 What is the format of the Fourier coefficients file used in Fourier coefficients file path: ? I don't see an example for this in svFSI-Tests. Can you send me an example file that I can use to test?

[vvedula22]

OK I will. You also asked for another example during the last solver meeting. Could you remind me what was that test case we talked about then? Will upload both these cases.

[CZHU20]

Hi @ktbolt, I saw you are creating templates for 2D/3D arrays. Is it possible to take advantage of existing libraries such as eigen3, boost etc to handle the matrix/vector operations for us?

[ktbolt]

@vvedula22 During the solver meeting I asked for a test case for variable wall properties.

[ktbolt]

@CZHU20 I wrote my own array and vector templates so I could easily reproduce Fortran arrays and debug indexing problems. I've been looking at Eigen and other packages for matrix/vector operations, would like to use something for the real svFSI implementation. I will probably use Eigen for the Python svZeroDSolver that I will help reimplement in C++.

What is your experience with matrix/vector packages? Do you have a preference?

[CZHU20]

Hi @ktbolt thanks for the explanation. I am only familiar with Eigen. Since it's just a bunch of header files, I think Eigen can also be easily packaged into svFSI for distribution.

[vvedula22]

@ktbolt Added two cases in svFSI-Tests :

• Fourier coefficients example is included in the pipe3D_RCR example
• CMM with variable wall properties is included in the vein-graft example

Also added README files to describe the problem for both these cases. Let me know if you have any questions on these.

[ktbolt]

Thanks @vvedula22 !

[ktbolt]

The way svFSI currently preprocesses input mesh, face, and BC data for a simulation is not very efficient, especially for an HPC cluster

• submit a job to the cluster
• wait for it to start
• execute svFSI using MPI
• process the mesh
• partition the mesh
• detect an error in the solver input file or other data
• job fails

Maybe what we need to do is to preprocess the input data and write that to a file that is guaranteed to be valid. This file can then be used for the simulation.

[ktbolt]

I discovered that I skipped converting the code that processes boundary conditions. This is probably important so I converted the 500 line READBC subroutine, added more XML parameters for BCs, code to read in various .dat files, etc.

[ktbolt]

Added the Add BoyFriend (Add_BF) XML parameter and converted the Fortran code processing this data.

[ktbolt]

I've finished converting the 400 line DISTRIBUTE function and its dependences that partitions meshes and distributes data to MPI processes.

So far I've converted 14 solver .inp files to XML, would like to convert the remaining 54 and test, make sure everything is working up to this stage in the conversion, always finding new parameters and a bug now and then.

[CZHU20]

@ktbolt Great achievement IMO! I have attempted a couple times, but never could make sense of how the jobs are split. Were you able to understand that fully?

[ktbolt]

Thanks @CZHU20 ! Converting SPLITJOBS took a lot of time to get right, no comments and the use of meaningless two-letter variables names. I understand what it is doing though, have written this sort of thing myself.

The use a key variable named j and comments like The work that suppose to be done by "l" reminds me of Men in Black ...

[ktbolt]

I've converted all of the solver input .inp files to .xml except for some of the higher-order elements tests in 02-lelas/01-beam_Euler-Bernoulli, found a few bugs and added more parameters and element types.

The 07-fsi/cmm/01-pipe_RCR/03-deform-wall_prestress/02-CMM-prestress test fails on MacOS 10.13 with a memory corruption error in the Fortran code using the latest svFSI release. It does not fail on Ubuntu 18.

The 09-shells/01-plate test produces VTK warnings

Building VTK_TRIANGLE 0 with less than three points, but VTK_TRIANGLE needs at least three points. Check the input.

in the C++ code when trying to read any of the face VTK files (.e.g. mesh/mesh-surfaces/x0.vtp). This is a 2D problem so the faces should be lines I think.

[CZHU20]

@ktbolt Hi Dave, I found a bug in the input file of 07-fsi/ale/02-channel-block-flag_2D. The fluid domain inlet was not fixed so it moved with the mesh motion. Just submitted a pull request, but I don't think I have the access to accept the request. You can see the file change here. https://github.com/SimVascular/svFSI-Tests/pull/5/commits

[vvedula22]

@ktbolt You may leave out the shell examples from testing for the time being until issue #46 is resolved. There could be some inconsistent output.

On the CMM case, I tried it on my Mac but could not reproduce the issue. However, I found a bug in the CMM code during memory allocation that might explain the error. Can you try applying the attached patch file and see if it fixes it? patch_bug_CMM_memalloc.txt

@CZHU20 I merged your pull request. Thanks for finding the bug.

[ktbolt]

@vvedula22 The 07-fsi/cmm/01-pipe_RCR/03-deform-wall_prestress/02-CMM-prestress test now works with your CMM fix. Thanks for the fix!

[ktbolt]

I am now converting INITIALIZE.f and the svFSILS code.

[ktbolt]

Found and fixed a few C++ compilation errors on Ubuntu.

[ktbolt]

I've fixed a bunch of array indexing problems in FSILS_LHS_CREATE. This guy seems to create data structures for manipulating the LHS matrix although it is not clear what the data structures are used for.

[CZHU20]

Hi Dave, I spent quite some time figuring out some of the data structures while implementing Hypre in svFSI. I am not sure what data structures you are referring to here, and as you probably have seen in the code, most of the data structures created in FSILS_LHS_CREATE are for parallel communication. Basically, there is the usual global-local indexing mapping and a local reordering that separates mesh points belonging solely to the current process with those shared with other processes. I have sent you a personal note. Hope it helps.

[ktbolt]

Thanks @CZHU20 for the information, very helpful! It is good to understand the general idea about what the code does.

[ktbolt]

@vvedula22 @CZHU20 I don't see any tests using the Simulation initialization file path keyword. Is this functionality that is useful?

[ktbolt]

I have added the Eigen C++ template library for linear algebra to svFSI/Code/ThirdParty. You can now use it like so

#include <eigen3/Eigen/Dense>

void compute_lu_invers()
{
  typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> MatrixType;

  MatrixType Ad(nd, nd);

  Eigen::FullPivLU<MatrixType> lu(Ad);

  lu.inverse();
}

However, I'm not using Eigen for any linear algebra operations, still using the Fortran LAPACK library, calling its Fortran functions from C++.

[ktbolt]

I've add a bunch of new code for material models, matrix operations, and post processing.

[ktbolt]

@vvedula22 I see that SUBROUTINE IFFT(gt, Y, dY) defined as

      PURE SUBROUTINE IFFT(gt, Y, dY)
      USE COMMOD

      IMPLICIT NONE
      TYPE(fcType), INTENT(IN) :: gt
      REAL(KIND=RKIND), INTENT(OUT) :: Y(gt%d), dY(gt%d)

is being called from SETBCDIRL using scalars rather than arrays like this

      REAL(KIND=RKIND) :: dirY, dirA

      ELSE IF (BTEST(lBc%bType,bType_ustd)) THEN
            CALL IFFT(lBc%gt, dirY, dirA)

I guess this might would work ok if gt%d is one. Is that the case?

[ktbolt]

I've finished converting the INITIALIZE function and all the functions called from that, added a bunch of functions from

txt.cpp
post.cpp
set_bc.cpp
output.cpp

I've not fully converted the code that writes text files, will do that later.

Found a few bugs and fixed them, problems with indexing into a subsection of an array, difficulty with variables named l and o and e ... and sometimes y.

The C++ code is about 35,000 lines right now.

[vvedula22]

@vvedula22 I see that SUBROUTINE IFFT(gt, Y, dY) defined as

      PURE SUBROUTINE IFFT(gt, Y, dY)
      USE COMMOD

      IMPLICIT NONE
      TYPE(fcType), INTENT(IN) :: gt
      REAL(KIND=RKIND), INTENT(OUT) :: Y(gt%d), dY(gt%d)

is being called from SETBCDIRL using scalars rather than arrays like this

      REAL(KIND=RKIND) :: dirY, dirA

      ELSE IF (BTEST(lBc%bType,bType_ustd)) THEN
            CALL IFFT(lBc%gt, dirY, dirA)

I guess this might would work ok if gt%d is one. Is that the case?

@ktbolt Do you still need my input on this one, or have you resolved this already?

You're right. In this particular function call, the dimension is 1 (gt%d=1) and so we can still make the function call. Fortran allows such calls as it always passes arguments reference. In other areas where the same function IFFT is called, you will see that an array is being passed. For e.g., see BF.f

[ktbolt]

@vvedula22 I've resolved this in the C++ conversion, was just curious about the Fortran implementation.

[ktbolt]

I've added support for Trilinos.

[CZHU20]

@vvedula22 @CZHU20 I don't see any tests using the Simulation initialization file path keyword. Is this functionality that is useful?

Hi @ktbolt Sorry for the late reply. I somehow missed the notification. Yes, this is a very useful functionality, especially in parametric study, in which I can use the results from one parameter as the initial condition for the other one to cut short the transient time to reach periodicity. The syntax would be Simulation initialization file path: ./result_1000.vtu. I have tested it and it worked fine in the current version.