Merge pull request #21 from keiyamamo/add_robinbc

add robin boundary condition

Author: dbruneau-mie

turtleFSI/modules/domain.py

@@ -4,10 +4,16 @@
 # PURPOSE.
 
 from turtleFSI.modules import *
-from dolfin import Constant, inner, inv, grad, div
+from dolfin import ds, MPI
 
+"""
+Date: 2022-11-07
+Comment from Kei: Naming here is a bit confusing. For example, dx_f, ds_s are hard to understand. 
+dx refers to the area of the domain while ds refers to the area of the boundary.
+f refere to the fluid while s refers to the solid.
+"""
 
-def assign_domain_properties(dx, dx_f_id, rho_f, mu_f, fluid_properties, dx_s_id, material_model, rho_s, mu_s, lambda_s, solid_properties, domains, **namespace):
+def assign_domain_properties(dx, dx_f_id, rho_f, mu_f, fluid_properties, dx_s_id, material_model, rho_s, mu_s, lambda_s, solid_properties, domains, ds_s_id, boundaries, robin_bc, **namespace):
     """
     Assigns solid and fluid properties to each region.   
     """
@@ -48,6 +54,19 @@ def assign_domain_properties(dx, dx_f_id, rho_f, mu_f, fluid_properties, dx_s_id
             solid_properties.append({"dx_s_id":dx_s_id,"material_model":material_model,"rho_s":rho_s,"mu_s":mu_s,"lambda_s":lambda_s})
     elif isinstance(solid_properties, dict): 
         solid_properties = [solid_properties]
+    
+    # RobinBC
+    if robin_bc:
+        ds_s = {}
+        if isinstance(ds_s_id, list): # If ds_s_id is a list (i.e, if there are multiple boundary regions):
+            for i, solid_boundaries in enumerate(ds_s_id):
+                ds_s[i] = ds(solid_boundaries, subdomain_data=boundaries) # Create ds_s for each boundary
+            ds_s_ext_id_list=ds_s_id
+        else:
+            ds_s[0] = ds(ds_s_id, subdomain_data=boundaries)
+            ds_s_ext_id_list=[ds_s_id] # If there aren't multpile boundary regions, and the boundary parameters are given as floats, convert solid parameters to lists.
+    else:
+        ds_s = None
+        ds_s_ext_id_list = None
 
-
-    return dict(dx_f=dx_f, dx_f_id_list=dx_f_id_list, fluid_properties=fluid_properties, dx_s=dx_s, dx_s_id_list=dx_s_id_list, solid_properties=solid_properties)
+    return dict(dx_f=dx_f, dx_f_id_list=dx_f_id_list, ds_s_ext_id_list=ds_s_ext_id_list, ds_s=ds_s, fluid_properties=fluid_properties, dx_s=dx_s, dx_s_id_list=dx_s_id_list, solid_properties=solid_properties)

turtleFSI/modules/solid.py

@@ -4,12 +4,11 @@
 # PURPOSE.
 
 from turtleFSI.modules import *
-
+from turtleFSI.problems import info_blue
 from dolfin import Constant, inner, grad, MPI
 
-
-def solid_setup(d_, v_, phi, psi, dx_s, dx_s_id_list, solid_properties, k, theta,
-                gravity, mesh, **namespace):
+def solid_setup(d_, v_, phi, psi, dx_s, ds_s, dx_s_id_list, ds_s_ext_id_list, solid_properties, k, theta,
+                gravity, mesh, robin_bc, k_s, c_s, **namespace):
 
     # DB added gravity in 3d functionality and multi material capability 16/3/21
     #
@@ -60,10 +59,23 @@ def solid_setup(d_, v_, phi, psi, dx_s, dx_s_id_list, solid_properties, k, theta
         # Stress (Note that if viscoelasticity is used, Piola1() is no longer the total stress, it is the non-rate dependant (elastic) component of the stress)
         F_solid_nonlinear += theta0 * inner(Piola1(d_["n"], solid_properties[solid_region]), grad(psi)) * dx_s[solid_region]
         F_solid_linear += theta1 * inner(Piola1(d_["n-1"], solid_properties[solid_region]), grad(psi)) * dx_s[solid_region]
-        # Gravity
+        # Gravity - y direction only
         if gravity is not None and mesh.geometry().dim() == 2:
             F_solid_linear -= inner(Constant((0, -gravity * rho_s)), psi)*dx_s[solid_region] 
         elif gravity is not None and mesh.geometry().dim() == 3:
             F_solid_linear -= inner(Constant((0, -gravity * rho_s,0)), psi)*dx_s[solid_region] 
+            
+    # Robin BC
+    """
+    The derivation comes from the eq.(9) in the followling paper:
+    Moireau, P., Xiao, N., Astorino, M. et al. External tissue support and fluid–structure simulation in blood flows. 
+    Biomech Model Mechanobiol 11, 1–18 (2012). https://doi.org/10.1007/s10237-011-0289-z
+    """
+    if robin_bc:
+        info_blue("Robin BC is used for the solid domain.")
+        for solid_boundaries in range(len(ds_s_ext_id_list)): 
+            F_solid_linear += theta0 * inner((k_s * d_["n"] + c_s * v_["n"]), psi)*ds_s[solid_boundaries] 
+            F_solid_linear += theta1 * inner((k_s * d_["n-1"] + c_s * v_["n-1"]), psi)*ds_s[solid_boundaries] 
+            
 
     return dict(F_solid_linear=F_solid_linear, F_solid_nonlinear=F_solid_nonlinear)

turtleFSI/problems/RobinBC_validation.py

@@ -0,0 +1,150 @@
+from dolfin import *
+from turtleFSI.problems import *
+import numpy as np
+from scipy.integrate import odeint
+from turtleFSI.utils.Probe import Probe
+import matplotlib.pyplot as plt
+
+"""
+This problem is a validation of the Robin BC implementation in the solid solver.
+The validation is done by using a mass-spring-damper system and comparing the results.
+We use cylinder mesh which is subjected to a constant gravity force in y-direction.
+This sciprt is meant to run with a single core, but can be easily parallelized.
+
+Mesh can be found in the following link:
+https://drive.google.com/drive/folders/1roV_iE_16Q847AQ_0tEsznIT-6EICX4o?usp=sharing
+"""
+
+# Set compiler arguments
+parameters["form_compiler"]["quadrature_degree"] = 6 
+parameters["form_compiler"]["optimize"] = True
+# The "ghost_mode" has to do with the assembly of form containing the facet normals n('+') within interior boundaries (dS). For 3D mesh the value should be "shared_vertex", for 2D mesh "shared_facet", the default value is "none".
+parameters["ghost_mode"] = "shared_vertex" #3D case
+_compiler_parameters = dict(parameters["form_compiler"])
+
+
+def set_problem_parameters(default_variables, **namespace):
+    # Overwrite default values
+    E_s_val = 1E6                              # Young modulus (elasticity) [Pa]
+    nu_s_val = 0.45                            # Poisson ratio (compressibility)
+    mu_s_val = E_s_val / (2 * (1 + nu_s_val))  # Shear modulus
+    lambda_s_val = nu_s_val * 2. * mu_s_val / (1. - 2. * nu_s_val)
+
+    # define and set problem variables values
+    default_variables.update(dict(
+        T=0.1,                              # Simulation end time
+        dt=0.0005,                          # Time step size
+        theta=0.501,                        # Theta scheme (implicit/explicit time stepping): 0.5 + dt
+        atol=1e-7,                          # Absolute tolerance in the Newton solver
+        rtol=1e-7,                          # Relative tolerance in the Newton solver
+        mesh_file="cylinder",               # Mesh file name
+        inlet_id=2,                         # inlet id 
+        outlet_id1=3,                       # outlet id
+        inlet_outlet_s_id=1011,             # solid inlet and outlet id
+        fsi_id=1022,                        # fsi Interface 
+        rigid_id=1011,                      # "rigid wall" id for the fluid and mesh problem
+        outer_wall_id=1033,                 # outer surface / external id 
+        ds_s_id=[1033],                     # ID of solid external wall (where we want to test Robin BC)
+        rho_f=1.025E3,                      # Fluid density [kg/m3]
+        mu_f=3.5E-3,                        # Fluid dynamic viscosity [Pa.s]
+        rho_s=1.0E3,                        # Solid density [kg/m3]
+        mu_s=mu_s_val,                      # Solid shear modulus or 2nd Lame Coef. [Pa]
+        nu_s=nu_s_val,                      # Solid Poisson ratio [-]
+        lambda_s=lambda_s_val,              # Solid 1rst Lamé coef. [Pa]
+        robin_bc = True,                    # Robin BC
+        k_s = 1.0E5,                        # elastic response necesary for RobinBC
+        c_s = 1.0E2,                        # viscoelastic response necesary for RobinBC 
+        dx_f_id=1,                          # ID of marker in the fluid domain
+        dx_s_id=1002,                       # ID of marker in the solid domain
+        extrapolation="laplace",            # laplace, elastic, biharmonic, no-extrapolation
+        extrapolation_sub_type="constant",  # constant, small_constant, volume, volume_change, bc1, bc2
+        recompute=30,                       # Number of iterations before recompute Jacobian. 
+        recompute_tstep=10,                 # Number of time steps before recompute Jacobian. 
+        save_step=1,                        # Save frequency of files for visualisation
+        folder="robinbc_validation",        # Folder where the results will be stored
+        checkpoint_step=50,                 # checkpoint frequency
+        kill_time=100000,                   # in seconds, after this time start dumping checkpoints every timestep
+        save_deg=1,                         # Default could be 1. 1 saves the nodal values only while 2 takes full advantage of the mide side nodes available in the P2 solution. P2 for nice visualisations 
+        gravity=2.0,                        # Gravitational force [m/s^2] 
+        fluid="no_fluid"                    # Do not solve for the fluid
+    ))
+
+    return default_variables
+
+
+def get_mesh_domain_and_boundaries(mesh_file, **namespace):
+    #Import mesh file
+    mesh = Mesh()
+    hdf = HDF5File(mesh.mpi_comm(), "mesh/" + mesh_file + ".h5", "r")
+    hdf.read(mesh, "/mesh", False)
+    boundaries = MeshFunction("size_t", mesh, 2)
+    hdf.read(boundaries, "/boundaries")
+    domains = MeshFunction("size_t", mesh, 3)
+    hdf.read(domains, "/domains")
+    
+    #Set all solid
+    domains.set_all(1002)
+
+    return mesh, domains, boundaries
+
+def create_bcs(**namespace):
+    """
+    In this problem we use Robin boundary condition which is implemented in the solid.py file.
+    Thus, we do not need specify any boundary conditions in this function. 
+    """
+    return dict(bcs=[])
+
+def _mass_spring_damper_system_ode(x, t, params_dict):
+    # Umpack parameters
+    F = params_dict['F'] # Volume of the domain   
+    A = params_dict['A'] # Area of the external surface (where Robin BC is applied)    
+    c = params_dict['c'] # Damping constant
+    k = params_dict['k'] # Stiffness of the spring 
+    m = params_dict['m'] # Mass of the domain
+    # Solve the system of ODEs
+    dx1dt = x[1]
+    dx2dt = (F - c*x[1]*A - k*x[0]*A)/m
+
+    dxdt = [dx1dt, dx2dt] 
+    return dxdt 
+
+def initiate(dvp_, mesh, DVP, **namespace):
+    # Position to probe
+    x_coordinate = mesh.coordinates()[:, 0]
+    y_coordinate = mesh.coordinates()[:, 1]
+    z_coordinate = mesh.coordinates()[:, 2]
+
+    x_middle = (x_coordinate.max() + x_coordinate.min())/2
+    y_middle = (y_coordinate.max() + y_coordinate.min())/2
+    z_middle = (z_coordinate.max() + z_coordinate.min())/2
+
+    middle_point = np.array([x_middle, y_middle, z_middle])
+    d_probe = Probe(middle_point, DVP.sub(0))
+    d_probe(dvp_["n"].sub(0, deepcopy=True))
+    return dict(d_probe=d_probe)
+
+def post_solve(d_probe, dvp_, **namespace):
+    d_probe(dvp_["n"].sub(0, deepcopy=True))
+
+def finished(T, dt, mesh, rho_s, k_s, c_s, boundaries, gravity, d_probe, **namespace):
+    # Define time step and initial conditions
+    t_analytical = np.linspace(0, T, int(T/dt) + 1)
+    analytical_solution_init = [0,0]
+    # Define parameters for the analytical solution
+    volume = assemble(1*dx(mesh))
+    ds_robin = Measure("ds", domain=mesh, subdomain_data=boundaries, subdomain_id=1033)
+    params_dict = dict()
+    params_dict["m"] = volume*rho_s
+    params_dict["k"] = k_s
+    params_dict["c"] = c_s
+    params_dict["A"] = assemble(1*ds_robin)
+    params_dict["F"] = -gravity*volume*rho_s
+    # Solve the ode to compute the analytical solution
+    analytical_solution = odeint(_mass_spring_damper_system_ode, analytical_solution_init, t_analytical, args=(params_dict,))
+    analytical_displacement = analytical_solution[:,0]
+    # Plot both numerical and analytical solutions for a comparison
+    plt.plot(d_probe.get_probe_sub(1), label="turtleFSI")
+    plt.plot(analytical_displacement, label="analytical")
+    plt.legend()
+    plt.show()
+