kvslab/turtlefsi

Monolithic Fluid-Structure Interaction (FSI) solver

turtlefsi's Issues

Strange numerical divergence in FSI example

Dear everyone,

To use water as the fluid medium, I made some changes to the benchmark case TF_fsi.
In order to keep the Reynolds number still be 200, I reduced the size of the entire computing domain by a factor of 100, and the inflow velocity $U_m$ by a factor of 10. This case can run well when the plate is rigid. However, when I added the solid plate part and activated the solid solver, some strange velocity and pressure surge emerged at the FSI surface near the free and fixed end of the plate, as you can see in the following figures. Do you know why this happened and how we can ovoid this problem?

The case and mesh file are also attached.

cylinderflap_tail_head_cy2_4D_bl.zip

Figure 1 Geometry and mesh

Figure 2 Results of rigid plate

Figure 3 Results of elastic plate

from dolfin import *
import numpy as np
from os import path
from mpi4py import MPI as pyMPI

from turtleFSI.problems import *
from turtleFSI.modules import *

parameters["ghost_mode"] = "shared_facet"
#_compiler_parameters = dict(parameters["form_compiler"])


def set_problem_parameters(default_variables, **namespace):
    # Overwrite or add new variables to 'default_variables'
    default_variables.update(dict(
        # Temporal variables
        T=2,                         # End time [s]
        dt=0.001,                      # Time step [s]
        theta=0.501,                    # Temporal scheme

        # Physical constants ('FSI 3')
        fluid_properties=[],
        solid_properties=[], 
        material_model="StVenantKirchoff",
        Um=0.2,                       # Max. velocity inlet, CDF3: 2.0 [m/s]
        rho_f=1.0e3,                  # Fluid density [kg/m3]
        mu_f=1.0e-3,                     # Fluid dynamic viscosity [Pa.s]
        rho_s=1.0e3,                  # Solid density[kg/m3]
        nu_s=0.4,                     # Solid Poisson ratio [-]
        mu_s=2.0e6,                   # Shear modulus, CSM3: 0.5E6 [Pa]
        lambda_s=4e6,                 # Solid 1st Lame Coefficient [Pa]
        fluid="fluid",                             # ["fluid", "no_fluid"] Turn off fluid and only solve the solid problem
        solid="solid",                             # ["solid", "no_solid"] Turn off solid and only solve the fluid problem


        # Problem specific
        folder="CylinderFlap-n2-tail-head-withCy-bl",      # Name of the results folder
        extrapolation="biharmonic",   # No displacement to extrapolate
        extrapolation_sub_type="constrained_disp",  # Biharmonic type
        bc_ids=[11, 12, 13, 16],          # Ids of makers for the mesh extrapolation
        #bc_ids=[2, 3, 4, 6],          # Ids of makers for the mesh extrapolation

        # Domain settings
        dx_f_id=17,       # Domain id of the fluid domain
        dx_s_id=18,       # Domain id of the solid domain

        # Solver settings
        recompute=1,                  # Compute the Jacobian matrix every iteration
        save_step=1,  # Save file frequency

        # Geometric variables
        R=0.0005,                       # Radius of the circle
        H=0.0041,                       # Total height
        L=0.25,                        # Length of domain
        f_L=0.0035,                     # Length of the flag
        f_H=0.0002,                     # Height of the flag
        c_x=0.002,                      # Center of the circle x-direction
        c_y=0.002))                     # Center of the circle y-direction

    default_variables["compiler_parameters"].update({"quadrature_degree": 5})

    return default_variables
"""
1 11 "Inlet"
1 12 "Outlet"
1 13 "Wall"
1 14 "Bar"
1 15 "Barwall"
1 16 "Circle"
2 17 "Fluid"
2 18 "Solid"
"""

def get_mesh_domain_and_boundaries(R, H, L, f_L, f_H, c_x, c_y, **namespace):
    # Read mesh
    mesh_folder = path.join(path.dirname(path.abspath(__file__)), "mesh")
    mesh_path = path.join(mesh_folder, "CylinderFlap-n2-tail-head-withCy2-4D-bl.h5")     # mesh geometry

    mesh = Mesh()

    hdf = HDF5File(mesh.mpi_comm(), mesh_path, "r")
    hdf.read(mesh, "/mesh", False)
    domains = MeshFunction("size_t", mesh, mesh.geometry().dim())
    hdf.read(domains, "/domains")
    boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)
    hdf.read(boundaries, "/boundaries")


    return mesh, domains, boundaries


def initiate(c_x, c_y, R, f_L, **namespace):
    # Coordinate for sampling statistics
    coord = [c_x + R + f_L+8*R, c_y] # point at the bar right end

    # Lists to hold results
    displacement_x_list = []
    displacement_y_list = []
    drag_list = []
    lift_list = []
    time_list = []

    return dict(displacement_x_list=displacement_x_list, displacement_y_list=displacement_y_list,
                drag_list=drag_list, lift_list=lift_list, time_list=time_list, coord=coord)


class Inlet(UserExpression):
    def __init__(self, Um, H, **kwargs):
        self.Um = 1.5 * Um
        self.H = H
        self.factor = 0
        super().__init__(**kwargs)

    def update(self, t):
        if t < 0.2:
            self.factor = 0.5 * (1 - np.cos(t * np.pi / 0.2)) * self.Um
        else:
            self.factor = self.Um

    def eval(self, value, x):
        value[0] = self.factor * x[1] * (self.H - x[1]) / (self.H / 2.0)**2
        value[1] = 0

    def value_shape(self):
        return (2,)


def create_bcs(DVP, v_deg, Um, H, boundaries, extrapolation_sub_type, **namespace):
    inlet = Inlet(Um, H, degree=v_deg)

    # Fluid velocity conditions
    u_inlet = DirichletBC(DVP.sub(1), inlet, boundaries, 11)
    u_wall = DirichletBC(DVP.sub(1), ((0.0, 0.0)), boundaries, 13)
    u_circ = DirichletBC(DVP.sub(1), ((0.0, 0.0)), boundaries, 16)
    u_barwall = DirichletBC(DVP.sub(1), ((0.0, 0.0)), boundaries, 15)

    # Pressure Conditions
    p_out = DirichletBC(DVP.sub(2), 0, boundaries, 12)

    # Assemble boundary conditions
    bcs = [u_wall, u_inlet, u_circ, u_barwall, p_out]

    # Boundary conditions on the displacement / extrapolation
    if extrapolation_sub_type != "constrained_disp_vel":
        d_wall = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 13)
        d_inlet = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 11)
        d_outlet = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 12)
        d_circle = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 16)
        d_barwall = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 15)
        for i in [d_wall, d_inlet, d_outlet, d_circle, d_barwall]:
            bcs.append(i)

    else:
        w_wall = DirichletBC(DVP.sub(3), ((0.0, 0.0)), boundaries, 13)
        w_inlet = DirichletBC(DVP.sub(3), ((0.0, 0.0)), boundaries, 11)
        w_outlet = DirichletBC(DVP.sub(3), ((0.0, 0.0)), boundaries, 12)
        w_circle = DirichletBC(DVP.sub(3), ((0.0, 0.0)), boundaries, 16)
        w_barwall = DirichletBC(DVP.sub(3), ((0.0, 0.0)), boundaries, 15)

        d_wall = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 13)
        d_inlet = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 11)
        d_outlet = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 12)
        d_circle = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 16)
        d_barwall = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 15)

        for i in [w_wall, w_inlet, w_outlet, w_circle, w_barwall,
                  d_wall, d_inlet, d_outlet, d_circle, d_barwall]:
            bcs.append(i)

    return dict(bcs=bcs, inlet=inlet)


def pre_solve(t, inlet, **namespace):
    """Update boundary conditions"""
    inlet.update(t)


################################################################################
# the function mpi4py_comm and peval are used to overcome FEniCS limitation of
# evaluating functions at a given mesh point in parallel.
# https://fenicsproject.discourse.group/t/problem-with-evaluation-at-a-point-in
# -parallel/1188


def mpi4py_comm(comm):
    '''Get mpi4py communicator'''
    try:
        return comm.tompi4py()
    except AttributeError:
        return comm


def peval(f, x):
    '''Parallel synced eval'''
    try:
        yloc = f(x)
    except RuntimeError:
        yloc = np.inf*np.ones(f.value_shape())

    comm = mpi4py_comm(f.function_space().mesh().mpi_comm())
    yglob = np.zeros_like(yloc)
    comm.Allreduce(yloc, yglob, op=pyMPI.MIN)

    return yglob
################################################################################


def post_solve(t, dvp_, coord, displacement_x_list, displacement_y_list, drag_list, lift_list, mu_f, n,
               verbose, time_list, ds, dS, **namespace):
    d = dvp_["n"].sub(0, deepcopy=True)
    v = dvp_["n"].sub(1, deepcopy=True)
    p = dvp_["n"].sub(2, deepcopy=True)

    # Compute drag and lift
    Dr = -assemble((sigma(v, p, d, mu_f)*n)[0]*ds(16))
    Li = -assemble((sigma(v, p, d, mu_f)*n)[1]*ds(16))
    Dr += -assemble((sigma(v("+"), p("+"), d("+"), mu_f)*n("+"))[0]*dS(14))
    Li += -assemble((sigma(v("+"), p("+"), d("+"), mu_f)*n("+"))[1]*dS(14))

    # Append results
    drag_list.append(Dr)
    lift_list.append(Li)
    time_list.append(t)
    d_eval = peval(d, coord)
    displacement_x_list.append(d_eval[0])
    displacement_y_list.append(d_eval[1])
    # displacement_x_list.append(d(coord)[0])
    # displacement_y_list.append(d(coord)[1])

    # Print
    if MPI.rank(MPI.comm_world) == 0 and verbose:
        print("Distance x: {:e}".format(displacement_x_list[-1]))
        print("Distance y: {:e}".format(displacement_y_list[-1]))
        print("Drag: {:e}".format(drag_list[-1]))
        print("Lift: {:e}".format(lift_list[-1]))


def finished(results_folder, displacement_x_list, displacement_y_list, drag_list, lift_list, time_list, **namespace):
    if MPI.rank(MPI.comm_world) == 0:
        np.savetxt(path.join(results_folder, 'Lift.txt'), lift_list, delimiter=',')
        np.savetxt(path.join(results_folder, 'Drag.txt'), drag_list, delimiter=',')
        np.savetxt(path.join(results_folder, 'Time.txt'), time_list, delimiter=',')
        np.savetxt(path.join(results_folder, 'dis_x.txt'), displacement_x_list, delimiter=',')
        np.savetxt(path.join(results_folder, 'dis_y.txt'), displacement_y_list, delimiter=',')

Error in make_womersley_bcs due to numpy.complex removal

When calling the make_womersley_bcs function in the Womersley.py module, I get the following error message:

AttributeError: module 'numpy' has no attribute 'complex'.

np.complex was a deprecated alias for the builtin complex and has been removed in recent versions of NumPy (see https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations). The suggested solution is to use complex instead of np.complex in the code. However, if the intention was to use the numpy scalar type, np.complex128 should be used instead. Which one should we use?

Can I simulate the self-propelled swimmer

Hi,

Thanks for your great tool.

Can I use the turtleFSI to simulate the self-propelled swimmer? How can I define the swimmer's behavior?

For example, I can prescribe the behavior by using a function of the midline of the body

Or specifically, can the turtleFSI be used to simulate the environment like in the paper “Reinforcement Learning and Wavelet Adapted Vortex Methods for Simulations of Self-propelled Swimmers”?

Thank you very much for your help.

Which boundaries should be chosen to be included in bc_ids[]?

Good day dear evryone:

I'm now trying to construct my own fsi case through gmsh to generate msh file to hdf5 file.
And that would need to use my own boundaries id number to set the case file.
Now I'am confused what boundaries should be chosen to be included in bc_ids[]?

In the TF_FSI case file, I think the [2,3,4,6] refer to inlet outlet wall and cylinder. Is my guess right?

# TF_FSI
        # Problem specific
        folder="TF_fsi_results",      # Name of the results folder
        extrapolation="biharmonic",   # No displacement to extrapolate
        extrapolation_sub_type="constrained_disp_vel",  # Biharmonic type
        bc_ids=[2, 3, 4, 6],          # Ids of makers for the mesh extrapolation
    
    # Mark boundaries
    Allboundaries = DomainBoundary() 
    boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)
    boundaries.set_all(0)        
    Allboundaries.mark(boundaries, 1)   
    Wall.mark(boundaries, 2)         
    Inlet.mark(boundaries, 3)       
    Outlet.mark(boundaries, 4)       
    Bar.mark(boundaries, 5)             
    Circle.mark(boundaries, 6)          
    Barwall.mark(boundaries, 7)         

and in the turtle_demo.py case file, it didn't define any bc_ids. that's confusing, because there are not the extrapolation_sub_type named constrained_disp in module file biharmonic.py, how can it know how to proceed with constrained_disp type?

 # turtle_demo.py
        extrapolation="biharmonic",    # Laplace, elastic, biharmonic, no-extrapolation
        extrapolation_sub_type="constrained_disp",  # ["constant", "small_constant", "volume", "volume_change", "constrained_disp", "constrained_disp_vel"]

 # biharmonic.py
    alpha_u = 0.01
    F_ext1 = alpha_u * inner(w_["n"], beta) * dx_f - alpha_u * inner(grad(d_["n"]), grad(beta)) * dx_f
    F_ext2 = alpha_u * inner(grad(w_["n"]), grad(phi)) * dx_f

    if extrapolation_sub_type == "constrained_disp_vel":
        for i in bc_ids:
            F_ext1 += alpha_u*inner(grad(d_["n"]) * n, beta) * ds(i)
            F_ext2 -= alpha_u*inner(grad(w_["n"]) * n, phi) * ds(i)

    F_fluid_linear += F_ext1 + F_ext2

    return dict(F_fluid_linear=F_fluid_linear)

Update documentation

Related to #47, we need to update the documentation, particularly related to installation.

$ conda create -n your_environment -c conda-forge turtleFSI

should be

conda env update --file environment.yml --name your_environment_name

and we should not use setup.py for building.

We should probably also consider moving from readdocs to github.io.

Fix drag and lift computations in TF problem

Interior facet integrals needs a volume tag to have a consistent orientation:

turtleFSI/turtleFSI/problems/TF_fsi.py

Lines 230 to 233 in e1b5350

	Dr = -assemble((sigma(v, p, d, mu_f)*n)[0]*ds(6))
	Li = -assemble((sigma(v, p, d, mu_f)*n)[1]*ds(6))
	Dr += -assemble((sigma(v("+"), p("+"), d("+"), mu_f)*n("+"))[0]*dS(5))
	Li += -assemble((sigma(v("+"), p("+"), d("+"), mu_f)*n("+"))[1]*dS(5))

ref
https://fenicsproject.discourse.group/t/integrating-over-an-interior-surface/247/4?u=dokken

Confirmation of a parameter of a standard case FSI3

Dear everyone,

I just found that there may be something confused me for the material settings of the benchmark case FSI3 in the case file TF_fsi.py
The first lame coefficient is set to be 4e6 but it may should be 8e6 from the paper.

        # Physical constants ('FSI 3')
        Um=2.0,                       # Max. velocity inlet, CDF3: 2.0 [m/s]
        rho_f=1.0e3,                  # Fluid density [kg/m3]
        mu_f=1.0,                     # Fluid dynamic viscosity [Pa.s]
        rho_s=1.0e3,                  # Solid density[kg/m3]
        nu_s=0.4,                     # Solid Poisson ratio [-]
        mu_s=2.0e6,                   # Shear modulus, CSM3: 0.5E6 [Pa]
        lambda_s=4e6,                 # Solid 1st Lame Coefficient [Pa]

Although the differences between 4e6 and 8e6 of first lame coefficient to calculate Elastic modulus is not too big, which are 5.333e6 and 5.6e6 pa.

Best regards
Zhang

mesh_path is overwritten

In PR #81, mesh_path is used as variable name in monolithic.py. This is a commonly used variable name in problem files and this change overwrites this variable resulting in errors when using the variable later in the problem file. Changing the variable name in monolithic.py to something else than mesh_path will fix this issue.

time loop conditional in `monolithic.py`

Inside monolithic.py, we have a conditional for time loop as follows.

turtleFSI/turtleFSI/monolithic.py

Line 156 in 0b95d44

while t <= T + dt / 10 and not stop: # + dt / 10 is a hack to ensure that we take the final time step t == T

I think + dt/10 should really be -dt/10 to ensure t==T because adding dt/10 makes it go into while loop even when t==T and this will end up t=T+dt as the final time step.

I wanna make sure that others agree on this, as this is important part of turtle.
Could I have your opinion on this?
@johannesring @jorgensd @dbruneau-mie

Creating a mesh with markers

Im having trouble creating a mesh

'get_mesh_domain_and_boundaries
The function is essential and unique for every problem, and has to be provided for the problem file to run. Here you read or define your mesh, domain markers, and boundary markers. In turtle_demo.py we read the mesh data from pre-existing “.xdmf” mesh files. In contrast to defining the domain and boundary markers using FEniCS functions, like in the the ‘Turek flag’-example (TF_fsi.py).

pygmsh and meshio are relevant tools to create geometries. For any question regarding meshing, we refer to the FEniCS documentation and discourse group.'

I want to create a simple circle and to visualise and calculate the lift drag and etc when air passes it

Some questions in the solid.py

Dear everyone

While I was reading the modules/solid.py , I got a few questions:

1.
The penalty function method is used to make sure v = d(d)/dt
And v = d(d)/dt has been changed to

This is really like the idea of the central difference.
And here, I don't know why to multiply phi(related to displacement) instead of psi(related to velocity)?
And as there has already a huge delta(1e7), will it have a difference if the parameter rho_s is removed?
Also this part represents rho_s*v_*δΔd is not fitting to the first Temporal term rho_s*a_*δΔv, how can the penalty part directly be added to the F_solid_linear?

        ## Temporal term and convection 
        F_solid_linear +=  (rho_s/k * inner(v_["n"] - v_["n-1"], psi)*dx_s[solid_region]
                                   + delta * rho_s * (1 / k) * inner(d_["n"] - d_["n-1"], phi)                * dx_s[solid_region]
                                    - delta * rho_s * inner(theta0 * v_["n"] + theta1 * v_["n-1"], phi) * dx_s[solid_region]) 

2.
for all the F linear part and nonlinear part in solid domain
The code includes the following:

        ## Temporal term and convection 
        F_solid_linear += (rho_s/k * inner(v_["n"] - v_["n-1"], psi)*dx_s[solid_region]
                          + delta * rho_s * (1 / k) * inner(d_["n"] - d_["n-1"], phi) * dx_s[solid_region]
                          - delta * rho_s * inner(theta0 * v_["n"] + theta1 * v_["n-1"], phi) * dx_s[solid_region]) # Maybe we could add viscoelasticity with v["n"] term instead of (1 / k) * inner(d_["n"] - d_["n-1"], phi) 
        
        # 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
        if gravity is not None and mesh.geometry().dim() == 2:
            F_solid_linear -= inner(Constant((0, -gravity * rho_s)), psi)*dx_s[solid_region] 

And I write it into formulas, and I also don't know what's the purpose of the term (1)(2)(3)(5) to multiply psi(related to velocity) instead of psi(related to velocity). I guess it is to use the concept of rate of work instead of work?

Any help would be greatly appreciated.

Best regards
ZHANG

Error in visualization of results in paraview (TF_fsi)

Hello !

I ran tutorial case TF_fsi with default parameters. I am facing issue while visualizing the velocty, pressure and displacement (.xdmf format) in paraview. I checked the data in another HDF5 viewer, it shows the data in array format. But, information of time is not perfectly visible. If possible, can I change writing file format ?

pytest missing laplace as mesh moving

I found out that the pytest actually does not test laplace.py because the parameter is not passed inside test_laplace function. The default parameter for extrapolation is set to biharmonic and I’m pretty sure that is the one used for test_laplace.

turtleFSI/tests/test_turtleFSI.py

Lines 67 to 86 in 29e8724

	@pytest.mark.parametrize("extrapolation_sub_type", ["volume", "volume_change",
	"constant", "small_constant"])
	def test_laplace(extrapolation_sub_type):
	cmd = ("turtleFSI --problem TF_fsi -dt 0.01 -T 0.05 --verbose True --theta 0.51" +
	" --folder tmp --sub-folder 4")
	d = system(cmd)
	drag = np.loadtxt("tmp/4/Drag.txt")[-1]
	lift = np.loadtxt("tmp/4/Lift.txt")[-1]
	distance_x = np.loadtxt("tmp/4/dis_x.txt")[-1]
	distance_y = np.loadtxt("tmp/4/dis_y.txt")[-1]
	distance_x_reference = -6.896013956339182e-06
	distance_y_reference = 1.876355330341896e-09
	drag_reference = 4.407481239804155
	lift_reference = -0.005404703556977697
	assert compare(distance_x, distance_x_reference)
	assert compare(distance_y, distance_y_reference)
	assert compare(drag, drag_reference)
	assert compare(lift, lift_reference)

Here is coverage report that shows we are not hitting laplace.py

--------- coverage: platform darwin, python 3.10.10-final-0 ----------
Name                                    Stmts   Miss  Cover
-----------------------------------------------------------
turtleFSI/__init__.py                       1      0   100%
turtleFSI/modules/__init__.py               1      0   100%
turtleFSI/modules/biharmonic.py            14      0   100%
turtleFSI/modules/common.py               111     69    38%
turtleFSI/modules/domain.py                43     22    49%
turtleFSI/modules/elastic.py               12      0   100%
turtleFSI/modules/fluid.py                 20      0   100%
turtleFSI/modules/laplace.py               16     16     0%
turtleFSI/modules/newtonsolver.py          52      2    96%
turtleFSI/modules/no_extrapolation.py       2      0   100%
turtleFSI/modules/no_fluid.py               8      0   100%
turtleFSI/modules/no_solid.py               8      0   100%
turtleFSI/modules/solid.py                 36     13    64%
turtleFSI/monolithic.py                   128     19    85%
turtleFSI/run_turtle.py                     8      1    88%
turtleFSI/utils/Probe.py                   70     49    30%
turtleFSI/utils/Womersley.py              138    119    14%
turtleFSI/utils/__init__.py                16     12    25%
turtleFSI/utils/argpar.py                 142     58    59%
-----------------------------------------------------------
TOTAL                                     826    380    54%

I will fix it and will create a PR.

Suspended solid in a rotating fluid - dealing with finite rotations

Physical scenario:
Consider a chunk of solid immersed in a fluid. The external fluid flow induces translations, rotations, and deformations on the solid.
As I understand it, in such a problem, the boundary conditions are imposed on the external fluid (for both velocity and pressure), and the displacement field solely relies on the extrapolation and uplifting methods used.

Let's focus on a rotating fluid, i.e., the boundary velocity is only forced to be in the angular direction. This fluid flow induces rotations of the solid. This external flow should cause rigid body rotation without significant deformation for a very rigid solid.

When simulating such a scenario in turtleFSI, the solid does rotate. However, this rotation also induces large deformations in the reference frame, causing it to "shrink," and eventually, these cause the simulation to diverge (I suspect it has to do with the jacobian diverging in the elements that shrink the most).

Concretely, using ParaView, here is a snapshot of the entire initial displacement field (wireframe, without cropping out the solid circle):

and here is the last step before the crash, with the WrapByVector filter (color code is the displacement field):

By now, I have tried both different constitutive relations (SVK, NH, and SVKEnergy), as well as different extrapolation methods (Laplace + volume change, and biharmonic + constrained_disp), but this issue persists.

Questions:

1. Is the ALE framework capable of solving such a problem where there are large rotations?
2. If so, how can turtleFSI deal with such a scenario? What would be possible solutions to this issue?

---

For simplicity, here are some parts of the code I use:

def set_problem_parameters(default_variables, **namespace):
    default_variables.update(dict(
        atol = 1e-5,                                # Absolute tolerance for the Newton solver
        rtol = 1e-5,                                 # Relative tolerance for the Newton solver
        T = 50,                                   # End time [s] (set T to several seconds to simulate several swimming cycles)
        tramp = 2.0,                                # Time to ramp-up to max inlet velocity (from 0 to Um)
        dt = 0.05,                                    # Time step [s]
        theta = 0.55,                              # Theta value (0.5 + dt), shifted Crank-Nicolson scheme
        Omega = 1,                             # Max. velocity inlet [m/s]
        pout = 0,                                 # Outlet pressure [Pa]  
        dx_f_id = 1,                           # Define two markers for the fluid domain
        rho_f = 1,                         # Fluid density [kg/m3]
        mu_f = 1,                        # Fluid dynamic viscosity [Pa*s]
        L = l,                                      # Experimental window [m]
        R = r,                                      # RBC radius [m]
        material_model="StVenantKirchoffEnergy",          # ["StVenantKirchoff", "NeoHookean", "StVenantKirchoffEnergy"]
        rho_s = 1,                                # Solid density [kg/m3]
        mu_s = 1,                                  # Solid shear modulus or 2nd Lame Coef. [Pa]
        lambda_s = 100,                             # Solid 1st Lame Coef. [Pa]
        nu_s = 1,                                  # Solid Poisson ratio [-]
        dx_s_id = 2,                                # ID of marker in the solid domain
        
        # Variational formulations
        fluid="fluid",                              # ["fluid", "no_fluid"] Turn off fluid and only solve the solid problem
        solid="solid",                              # ["solid", "no_solid"] Turn off solid and only solve the fluid problem
        extrapolation = "biharmonic",                  # Laplace, elastic, biharmonic, no-extrapolation
        extrapolation_sub_type = "constrained_disp",   # ["constant", "small_constant", "volume", "volume_change", "constrained_disp", "constrained_disp_vel"]
        
        recompute = 20,                              # Recompute the Jacobian matrix every "recompute" Newton iterations
        N = 21,                                     # number of points along x or y axis when creating structured mesh
        checkpoint_step = 10,                       # Store results for restart every 10th timestep
        folder = "Rotation_results"),              # Name of the folder to save the data
        save_step = 5,                                 # Frequency of data saving
        
        # Spatial settings
        lmbda = 1.0,                                # Relaxation parameter for the extrapolation
        max_it = 100,                               # Maximum number of Newton iterations
        v_deg = 2,                                    # velocity degree
        p_deg = 1,                                    # pressure degree
        d_deg = 2                                     # solid deformation degree
    )
    default_variables["compiler_parameters"].update({"quadrature_degree": 5})

    return default_variables

Here is my mesh construction (a fluid square, with a circular solid inside):

def get_mesh_domain_and_boundaries(L, R,N, **namespace):
    
    # # Mesh and boundaries for square domain
    mesh = RectangleMesh(Point(-L, -L), Point(L,L), N, N, "crossed")
    BottomWall = AutoSubDomain(lambda x: near(x[1], -L))
    TopWall = AutoSubDomain(lambda x: near(x[1], L))
    LeftWall = AutoSubDomain(lambda x: near(x[0], -L))
    RightWall = AutoSubDomain(lambda x: near(x[0], L))
    
    def capsule(x):
        return ((x[0])**2 + (x[1] )**2 < (((R)**2)))

    # Mark boundaries
    Allboundaries = DomainBoundary()
    boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)
    boundaries.set_all(0)
    Allboundaries.mark(boundaries, 1)
    
    # Mark square boundaries
    TopWall.mark(boundaries, 2)
    RightWall.mark(boundaries, 3)
    BottomWall.mark(boundaries, 4)
    LeftWall.mark(boundaries, 5)

    # Define the domain
    CapsuleArea = AutoSubDomain(capsule)
    
    domains = MeshFunction("size_t", mesh, mesh.geometry().dim())
    domains.set_all(1)
    CapsuleArea.mark(domains, 2)
    
    return mesh, domains, boundaries

I define the velocity at the boundary as

class ExtRotation(UserExpression):
    def __init__(self, Omega, tramp, L, **kwargs):
        self.t = 0.0
        self.t_ramp = tramp  # time to ramp-up to max inlet velocity (from 0 to Um)
        self.Omega = Omega       # Max. velocity inlet [m/s]
        self.L = L
        super().__init__(**kwargs)

    def update(self, t):
        self.t = t
        self.value = self.Omega

    def eval(self, value, x):
        value[0] = self.value * x[1] / self.L
        value[1] = -self.value * x[0] / self.L

    def value_shape(self):
        return (2,)

Then, I have for the boundary conditions

def create_bcs(DVP, boundaries, L,Omega, pout, tramp, v_deg, p_deg, extrapolation_sub_type, verbose, **namespace):
    if MPI.rank(MPI.comm_world) == 0 and verbose:
        print("Create bcs")
    
    extRotation = ExtRotation(Omega, tramp, L, degree=v_deg)
    
    # Square geometry boundary conditions
    # External shear fluid velocity conditions 
    u_left = DirichletBC(DVP.sub(1), extRotation, boundaries, 5)
    u_right = DirichletBC(DVP.sub(1), extRotation, boundaries, 3)
    u_top = DirichletBC(DVP.sub(1), extRotation, boundaries, 2)
    u_bot = DirichletBC(DVP.sub(1), extRotation, boundaries, 4)
    
    # Pressure boundary conditions
    p_left = DirichletBC(DVP.sub(2), pout, boundaries, 5)
    p_right = DirichletBC(DVP.sub(2), pout, boundaries, 3)
    p_top = DirichletBC(DVP.sub(2), pout, boundaries, 2)
    p_bot = DirichletBC(DVP.sub(2), pout, boundaries, 4)
    
    # Assemble boundary conditions
    bcs = [u_top, u_right, u_bot, u_left, p_left, p_bot, p_right, p_top]
    
    return dict(bcs=bcs, extRotation=extRotation, extRotationDsip = extRotationDsip)

And I also have to update the boundary velocity at each timestep using pre_solve:

def pre_solve(t, extRotation, **namespace): 
    extRotation.update(t)

3D FSI pipe flow problem

Dears,
I can't find the ' 3D FSI pipe flow problem' in the problem folder.
Please guide me to where that could be found.

Questions regarding Python versions

In the documentation it is stated that Python >=3.5 is required.
However, to use the full version of turtleFSI, fenicstools and cppimport is required.
cppimport does not support Python < 3.7.

Should the documentation really state that Python >=3.7 is required?

`setup.py` not working properly

When I run

python setup.py install

on my local machine or cluster, I get the following error message

No local packages or working download links found for fenics-dolfin
error: Could not find suitable distribution for Requirement.parse('fenics-dolfin')

I believe it is due to the recent update at the following line.

turtleFSI/setup.py

Line 6 in e1b5350

DEPENDENCIES = ['configargparse', "fenics-dolfin",

I remember it worked fine without this fenics-dolfin dependencies before, so I’m wondering if we need it.
Could @jorgensd comment on this?

Rotating a flap around a point on a aerofoil every timestamp

Hi,
now successfully made the mesh for aerofoil and positioned the flap on the point where it should hinge. Now i want the flap to rotate about the hinge for every timestamp. But I incorporated it but it remaining the same. My flap angle is changing but mesh isnt rotating.
i put it from -20 to 160 degreed as that when it will be parallel to the aerofoil. Below is the image when andle is -20. so when its 160 it would have rotated 180 degrees anticlockwise.

heres my code:

File under GNU GPL (v3) licence, see LICENSE file for details.

This software is distributed WITHOUT ANY WARRANTY; without even

the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR

PURPOSE.

"""
[1] Turek, Stefan, and Jaroslav Hron. "Proposal for numerical benchmarking of fluid-structure interaction
between an elastic object and laminar incompressible flow." Fluid-structure interaction.
Springer, Berlin, Heidelberg, 2006. 371-385."""
import mshr
from dolfin import *
import numpy as np
from os import path

from turtleFSI.problems import *
from turtleFSI.modules import *

def set_problem_parameters(default_variables, **namespace):
# Overwrite or add new variables to 'default_variables'
default_variables.update(dict(
# Temporal variables
T=0.5, # End time [s]
dt=0.01, # Time step [s]
theta=0.5, # Temporal scheme

    # Physical constants
    Um=1.0,                       # Max. velocity inlet, CDF3: 2.0 [m/s]
    rho_f=1.0e3,                  # Fluid density [kg/m3]
    mu_f=1.0,                     # Fluid dynamic viscosity [Pa.s]
    rho_s=1.0e3,                  # Solid density[kg/m3]
    nu_s=0.4,                     # Solid Poisson ratio [-]
    mu_s=2.0e6,                   # Shear modulus, CSM3: 0.5E6 [Pa]
    lambda_s=4e6,                 # Solid 1st Lame Coefficient [Pa]

    # Problem specific
    folder="ellipse_results",# Name of the results folder
    fluid="fluid",
    solid="no_solid",         # Do not solve for the solid
    extrapolation="no_extrapolation",  # No displacement to extrapolate
    current_flap_angle = -20,
    H=0.8,                        # Total height
    L=2.2,                        # Length of domain, adjusted to fit airfoil
    airfoil_chord_length=0.1,     # Chord length of the airfoil, adjust as needed
    airfoil_position=(0.6, 0.4),  # Position of the airfoil within the domain, adjust as needed
))

return default_variables

def rotate_point(point, angle_deg, rotation_center=(0, 0)):
"""Rotate a point counterclockwise by a given angle around a rotation center."""
angle_rad = np.radians(angle_deg)
ox, oy = rotation_center
px, py = point

qx = ox + np.cos(angle_rad) * (px - ox) - np.sin(angle_rad) * (py - oy)
qy = oy + np.sin(angle_rad) * (px - ox) + np.cos(angle_rad) * (py - oy)

return qx, qy

def rotate_shape(points, angle_deg, rotation_center):
"""Rotate a shape counterclockwise by a given angle around a rotation center."""
rotated_points = [rotate_point(pt, angle_deg, rotation_center) for pt in points]
return rotated_points

def generate_airfoil(chord_length, num_points=100, angle_of_attack=-20, translate_x=0.2, translate_y=0.5):
"""Generate NACA 0012 airfoil points and translate them."""
m, p, t = 0, 0, 0.12 # NACA 0012 parameters
x = np.linspace(0, 1, num_points)
yt = 5 * t * (0.2969 * np.sqrt(x) - 0.1260 * x - 0.3516 * x2 + 0.2843 * x3 - 0.1015 * x**4)

airfoil_points = [(xi, yi) for xi, yi in zip(x, yt)] + [(xi, -yi) for xi, yi in zip(reversed(x), reversed(yt))]

# Rotate airfoil
rotated_airfoil_points = [rotate_point(pt, angle_of_attack) for pt in airfoil_points]

# Translate airfoil
translated_airfoil_points = [(x + translate_x, y + translate_y) for x, y in rotated_airfoil_points]

return translated_airfoil_points

def sort_points_counter_clockwise(points):
if not points:
# Handle empty points array
return []

# Compute centroid of the points
points_array = np.array(points)
if points_array.ndim != 2 or points_array.shape[1] != 2:
    # Handle incorrect points structure
    raise ValueError("Points array is incorrectly structured.")

centroid = np.mean(points_array, axis=0)
cx, cy = centroid

# Sort points based on the angle from the centroid
sorted_points = sorted(points, key=lambda point: np.arctan2(point[1] - cy, point[0] - cx))

return sorted_points

def create_rectangular_flap(airfoil_points, flap_width, flap_length_ratio, position_ratio, vertical_offset):
chord_length = max(pt[0] for pt in airfoil_points) - min(pt[0] for pt in airfoil_points)
flap_length = chord_length * flap_length_ratio
flap_start_x = chord_length * position_ratio + min(pt[0] for pt in airfoil_points)

# Define the four corners of the flap rectangle
bottom_left = (flap_start_x, -flap_width / 2 + vertical_offset)
bottom_right = (flap_start_x + flap_length, -flap_width / 2 + vertical_offset)
top_right = (flap_start_x + flap_length, flap_width / 2 + vertical_offset)
top_left = (flap_start_x, flap_width / 2 + vertical_offset)

return [bottom_left, bottom_right, top_right, top_left]

def near_boundary(point, boundary_points, tolerance=0.005):
"""
Check if a point is near the given set of boundary points.
"""
return any(np.linalg.norm(np.array(point) - np.array(boundary_point)) < tolerance for boundary_point in boundary_points)

class AirfoilBoundary(SubDomain):
def init(self, sorted_main_airfoil_points, tolerance=0.005):
super().init()
self.sorted_main_airfoil_points = sorted_main_airfoil_points
self.tolerance = tolerance

def inside(self, x, on_boundary):
    return on_boundary and near_boundary(x, self.sorted_main_airfoil_points, self.tolerance)

class FlapBoundary(SubDomain):
def init(self, sorted_flap_points,sorted_main_airfoil_points, tolerance=0.005):
super().init()
self.sorted_flap_points = sorted_flap_points
self.sorted_main_airfoil_points = sorted_main_airfoil_points
self.tolerance = tolerance

def inside(self, x, on_boundary):
    is_on_flap = near_boundary(x, self.sorted_flap_points, self.tolerance)
    is_not_on_airfoil = not near_boundary(x, self.sorted_main_airfoil_points, self.tolerance)
    return on_boundary and is_on_flap and is_not_on_airfoil

**def update_flap_angle(current_angle, min_angle=-20, max_angle=160, step=1):
"""
Update the flap angle within a specified range.

Parameters:
- current_angle: The current angle of the flap.
- min_angle: The minimum angle limit.
- max_angle: The maximum angle limit.
- step: The step size for the angle update.

Returns:
- The updated angle.
"""
# Increment or decrement the angle based on its current value
if current_angle >= max_angle or current_angle < min_angle:
    step *= -1  # Change direction at extremes

new_angle = current_angle + step

# Ensure the new angle stays within the bounds
if new_angle > max_angle:
    new_angle = max_angle
elif new_angle < min_angle:
    new_angle = min_angle

return new_angle**

def get_mesh_domain_and_boundaries(L,H,current_flap_angle,**namespace):
channel = mshr.Rectangle(Point(0.0, 0.0), Point(L, H))
main_airfoil_points = generate_airfoil(chord_length=8, angle_of_attack=-20, translate_x=0.2, translate_y=0.5)
sorted_main_airfoil_points = sort_points_counter_clockwise(main_airfoil_points)
vertical_offset = 0.345
flap_width = 0.01
flap_points = create_rectangular_flap(main_airfoil_points, flap_width, 0.2, 0.6, vertical_offset)
rotation_center = min(flap_points, key=lambda pt: pt[0]) # Hinge point
rotated_flap_points = [rotate_point(pt, current_flap_angle, rotation_center) for pt in flap_points]
sorted_flap_points = sort_points_counter_clockwise(rotated_flap_points)
main_airfoil = mshr.Polygon([Point(x, y) for x, y in sorted_main_airfoil_points])
flap = mshr.Polygon([Point(x, y) for x, y in sorted_flap_points])

domain = channel - main_airfoil - flap
mesh = mshr.generate_mesh(domain, 100)


# Define the boundaries
inlet_boundary = AutoSubDomain(lambda x: near(x[0], 0))  # Inlet at x = 0
outlet_boundary = AutoSubDomain(lambda x: near(x[0], L))  # Outlet at x = L
wall_boundary = AutoSubDomain(lambda x: near(x[1], 0) or near(x[1], H))  # Walls at y = 0 and y = H
airfoil_boundary = AirfoilBoundary(sorted_main_airfoil_points)
flap_boundary = FlapBoundary(sorted_flap_points,sorted_main_airfoil_points)  

# Mark the boundaries
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
inlet_boundary.mark(boundaries, 1)
wall_boundary.mark(boundaries, 2)
outlet_boundary.mark(boundaries, 3)
airfoil_boundary.mark(boundaries, 4)  # Assuming 4 is the airfoil boundary
flap_boundary.mark(boundaries, 5)  # Assuming 5 is the flap boundary

# Define the domain
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(1)

return mesh, domains, boundaries

def initiate(**namespace):
# Lists to hold displacement, forces, and time
drag_list = []
lift_list = []
time_list = []

return dict(drag_list=drag_list, lift_list=lift_list, time_list=time_list)

class Inlet(UserExpression):
def init(self, Um, H, **kwargs):
self.Um = Um
self.H = H
self.factor = 0
super().init(**kwargs)

def update(self, t):
        self.Um =  self.Um

def eval(self, value, x):
    value[0] = self.Um  # Set the x-component of the velocity to the freestream velocity
    value[1] = 0  # y-component of the velocity is zero

def value_shape(self):
    return (2,)

def create_bcs(DVP, v_deg, Um, H, boundaries, **namespace):
inlet_velocity = Inlet(Um, H, degree=v_deg)
u_inlet = DirichletBC(DVP.sub(1), inlet_velocity, boundaries, 1)

u_wall = DirichletBC(DVP.sub(1), ((0.0, 0.0)), boundaries, 2)

u_airfoil = DirichletBC(DVP.sub(1), ((0.0, 0.0)), boundaries, 4)  # Airfoil boundary

u_flap = DirichletBC(DVP.sub(1), ((0.0, 0.0)), boundaries, 5)  # Flap boundary

p_out = DirichletBC(DVP.sub(2), 0, boundaries, 3)  # Outlet pressure

bcs = [u_inlet, u_wall, u_airfoil, u_flap, p_out]

return dict(bcs=bcs, inlet=inlet_velocity)

Example usage in the simulation loop

def pre_solve(t, inlet, **namespace):
inlet.update(t)

def post_solve(t, dvp_, n, L, H,DVP, v_deg, Um, boundaries,current_flap_angle, drag_list, lift_list, time_list, mu_f, verbose, ds, namespace):
# Get deformation, velocity, and pressure
d = dvp_["n"].sub(0, deepcopy=True)
v = dvp_["n"].sub(1, deepcopy=True)
p = dvp_["n"].sub(2, deepcopy=True)
# Compute forces
force = dot(sigma(v, p, d, mu_f), n)
drag_list.append(-assemble(force[0]*ds(4)))
lift_list.append(-assemble(force[1]*ds(4)))
time_list.append(t)
# Print results
if MPI.rank(MPI.comm_world) == 0 and verbose:
print("Drag: {:e}".format(drag_list[-1]))
print("Lift: {:e}".format(lift_list[-1]))
**current_flap_angle = update_flap_angle(current_flap_angle)
print(f"Current Flap Angle: {current_flap_angle}")
new_mesh, new_domains, new_boundaries = get_mesh_domain_and_boundaries(L, H, current_flap_angle)
return new_mesh, new_domains, new_boundaries, current_flap_angle **

def finished(drag_list, lift_list, time_list, results_folder, **namespace):
# Store results when the computation is finished
if MPI.rank(MPI.comm_world) == 0:
np.savetxt(path.join(results_folder, 'Lift.txt'), lift_list, delimiter=',')
np.savetxt(path.join(results_folder, 'Drag.txt'), drag_list, delimiter=',')
np.savetxt(path.join(results_folder, 'Time.txt'), time_list, delimiter=',')

and this is how i adapted monolithic:

while t <= T + dt / 10 and not stop: # + dt / 10 is a hack to ensure that we take the final time step t == T
t += dt

# Pre solve hook
tmp_dict = pre_solve(**vars())
if tmp_dict is not None:
    vars().update(tmp_dict)

# Solve
vars().update(newtonsolver(**vars()))

# Update vectors
for i, t_tmp in enumerate(times[:-1]):
    dvp_[t_tmp].vector().zero()
    dvp_[t_tmp].vector().axpy(1, dvp_[times[i+1]].vector())

# After solve hook
mesh, domains, boundaries, current_flap_angle = post_solve(**vars())
# Checkpoint
if counter % checkpoint_step == 0:
    checkpoint(**vars())

# Store results
if counter % save_step == 0:
    vars().update(save_files_visualization(**vars()))

Making a mesh

Hi, i was wondering how to make a mesh for my turtleFSI
I edited TF_cfd and removed the flag properties but need to edit the mesh to do the same but don't now how to adjust it .

File under GNU GPL (v3) licence, see LICENSE file for details.

This software is distributed WITHOUT ANY WARRANTY; without even

the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR

PURPOSE.

"""Problem file for running the "CFD" benchmarks in [1]. The problem is a channel flow
with a circle and a flag attached to it. For the CFD problem both the circle and flag is rigid.

[1] Turek, Stefan, and Jaroslav Hron. "Proposal for numerical benchmarking of fluid-structure interaction
between an elastic object and laminar incompressible flow." Fluid-structure interaction.
Springer, Berlin, Heidelberg, 2006. 371-385."""

from dolfin import *
import numpy as np
from os import path

from turtleFSI.problems import *
from turtleFSI.modules import *

def set_problem_parameters(default_variables, **namespace):
# Overwrite or add new variables to 'default_variables'
default_variables.update(dict(
# Temporal variables
T=30, # End time [s]
dt=0.01, # Time step [s]
theta=0.5, # Temporal scheme

    # Physical constants
    rho_f=1.0E3,              # Fluid density [kg/m3]
    mu_f=1.0,                 # Fluid dynamic viscosity [Pa.s]
    Um=2.0,                   # Max. velocity inlet (CDF3: 2.0) [m/s]

    # Problem specific
    folder="TF_cfd_results",  # Name of the results folder
    solid="no_solid",         # Do not solve for the solid
    extrapolation="no_extrapolation",  # No displacement to extrapolate

    # Geometric variables
    H=0.41,                   # Total height
    L=2.5))                   # Length of domain

return default_variables

def get_mesh_domain_and_boundaries(L, H, **namespace):
# Load and refine mesh
mesh = Mesh(path.join(path.dirname(path.abspath(file)), "..", "mesh", "TF_cfd.xml.gz"))
mesh = refine(mesh)

# Define the boundaries
Inlet = AutoSubDomain(lambda x: near(x[0], 0))
Outlet = AutoSubDomain(lambda x: near(x[0], L))
Walls = AutoSubDomain(lambda x: near(x[1], 0) or near(x[1], H))

# Mark the boundaries
Allboundaries = DomainBoundary()
boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)
boundaries set_all(0)
Allboundaries.mark(boundaries, 4)  # Circle and flag
Inlet.mark(boundaries, 1)
Walls.mark(boundaries, 2)
Outlet.mark(boundaries, 3)

# Define the domain
domains = MeshFunction("size_t", mesh, mesh.geometry().dim())
domains.set_all(1)

return mesh, domains, boundaries

def initiate(**namespace):
# Lists to hold displacement, forces, and time
drag_list = []
lift_list = []
time_list = []

return dict(drag_list=drag_list, lift_list=lift_list, time_list=time_list)

class Inlet(UserExpression):
def init(self, Um, H, **kwargs):
self.Um = Um * 1.5
self.H = H
self.factor = 0
super().init(**kwargs)

def update(self, t):
    if t < 2:
        self.factor = 0.5 * (1 - np.cos(t * np.pi / 2)) * self.Um
    else:
        self.factor = self.Um

def eval(self, value, x):
    value[0] = self.factor * x[1] * (self.H - x[1]) / (self.H / 2.0)**2
    value[1] = 0

def value_shape(self):
    return (2,)

def create_bcs(DVP, Um, H, v_deg, boundaries, **namespace):
# Create inlet expression
inlet = Inlet(Um, H, degree=v_deg)

# Fluid velocity conditions
u_inlet = DirichletBC(DVP.sub(1), inlet, boundaries, 1)
u_wall = DirichletBC(DVP.sub(1), ((0.0, 0.0)), boundaries, 2)

# Pressure Conditions
p_out = DirichletBC(DVP.sub(2), 0, boundaries, 3)

return dict(bcs=[u_wall, u_inlet, p_out], inlet=inlet)

def pre_solve(t, inlet, **namespace):
"""Update boundary conditions"""
inlet.update(t)

def post_solve(t, dvp_, n, drag_list, lift_list, time_list, mu_f, verbose, ds, **namespace):
# Get deformation, velocity, and pressure
d = dvp_["n"].sub(0, deepcopy=True)
v = dvp_["n"].sub(1, deepcopy=True)
p = dvp_["n"].sub(2, deepcopy=True)

# Compute forces
force = dot(sigma(v, p, d, mu_f), n)
drag_list.append(-assemble(force[0]*ds(4)))
lift_list.append(-assemble(force[1]*ds(4)))
time_list.append(t)

# Print results
if MPI.rank(MPI.comm_world) == 0 and verbose:
    print("Drag: {:e}".format(drag_list[-1]))
    print("Lift: {:e}".format(lift_list[-1]))

def finished(drag_list, lift_list, time_list, results_folder, **namespace):
# Store results when the computation is finished
if MPI.rank(MPI.comm_world) == 0:
np.savetxt(path.join(results_folder, 'Lift.txt'), lift_list, delimiter=',')
np.savetxt(path.join(results_folder, 'Drag.txt'), drag_list, delimiter=',')
np.savetxt(path.join(results_folder, 'Time.txt'), time_list, delimiter=',')

test_save_deg2 fails for FEniCS stable

The test test_save_deg2 fails for FEniCS stable with the following output:

Newton iteration 0: r (atol) = 4.276e+02 (tol = 1.000e-07), r (rel) = 5.250e+01 (tol = 1.000e-07)
Traceback (most recent call last):
  File "/usr/share/miniconda3/envs/turtleFSI/bin/turtleFSI", line 8, in <module>
    sys.exit(main())
  File "/usr/share/miniconda3/envs/turtleFSI/lib/python3.8/site-packages/turtleFSI/run_turtle.py", line 18, in main
    from turtleFSI import monolithic
  File "/usr/share/miniconda3/envs/turtleFSI/lib/python3.8/site-packages/turtleFSI/monolithic.py", line 195, in <module>
    vars().update(save_files_visualization(**vars()))
  File "/usr/share/miniconda3/envs/turtleFSI/lib/python3.8/site-packages/turtleFSI/problems/__init__.py", line 300, in save_files_visualization
    namespace["d_viz"].vector()[:] = namespace["dv_trans"]*d.vector()
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
Compute Jacobian matrix
*** using the information listed below, you can ask for help at
***
***     [email protected]
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'MatCreateVecs'.
*** Reason:  PETSc error code is: 98 (General MPI error).
*** Where:   This error was encountered inside /home/conda/feedstock_root/build_artifacts/fenics-pkgs_1696906530109/work/dolfin/dolfin/la/PETScBaseMatrix.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset:  2e001bd1aae8e14d758264f77382245e6eed04b0
*** -------------------------------------------------------------------------

Newton iteration 1: r (atol) = 5.941e-03 (tol = 1.000e-07), r (rel) = 9.318e+01 (tol = 1.000e-07)
Compute Jacobian matrix
Newton iteration 2: r (atol) = 2.873e-04 (tol = 1.000e-07), r (rel) = 1.181e-04 (tol = 1.000e-07)
Compute Jacobian matrix
Newton iteration 3: r (atol) = 2.515e-11 (tol = 1.000e-07), r (rel) = 3.317e-11 (tol = 1.000e-07)
Distance x: -2.527795e-06
Distance y: -3.826376e-10
Drag: 1.088597e+00
Lift: -9.584695e-04
FAILED

The test only fails for ubuntu-latest and for Python 3.8 and 3.10, while the test runs fine for Python 3.9 and 3.11.

The complete output is available here.

Definitions of some quantities inside newtonsolver

Dear everyone:

The core calculating part inside newtonsolver writes:

        if recompute_for_timestep or recompute_frequency or recompute_residual or recompute_initialize:
            if MPI.rank(MPI.comm_world) == 0 and verbose:
                print("Compute Jacobian matrix")
            A = assemble(J_nonlinear, tensor=A,
                         form_compiler_parameters=compiler_parameters,
                         keep_diagonal=True)
            A.axpy(1.0, A_pre, True)
            A.ident_zeros()
            [bc.apply(A) for bc in bcs]
            up_sol.set_operator(A)

        # Compute right hand side
        b = assemble(-F, tensor=b)

        # Apply boundary conditions and solve
        [bc.apply(b, dvp_["n"].vector()) for bc in bcs]
        up_sol.solve(dvp_res.vector(), b)
        dvp_["n"].vector().axpy(lmbda, dvp_res.vector())
        [bc.apply(dvp_["n"].vector()) for bc in bcs]

And in my understanding, for turtlefsi, in each time step, we will calculate the increments of displacement, velocity and pressure and then add the increments part to the old displacement, velocity and pressure.
Does it mean:
trialfunction chi mean real Δd Δv Δp
testfunction psi phi gamma mean virtual δ(Δd) δ(Δv) δ(Δp)
Most importantly, if dvp_res is for storing the results of real Δd Δv Δp (the increments of displacement, velocity and pressure)

Only one thing I am sure about:
dvp_["n"] means displacement, velocity and pressure at time step n,
dvp_["n-1"] means displacement, velocity and pressure at time step n-1,

Do I understand right for these quantities?

Best regards
Zhang

Issue new release of TurtleFSI to generate DOI for Nature Comms Paper

Hi @keiyamamo @jorgensd @johannesring @KristianValen-Sendstad, I would like to issue a new release (v2.2.1) of TurtleFSI to generate a DOI in Zenodo, so I can reference it in my Nature Comms paper. Do you have any concerns about doing this?

Remove parameterization over conda versions

Unfortunately, conda updates python version when we try to pin it, resulting in unnecessary usage of CI.

Try conda only with latest (python 3.10 or python 3.11)

Does this open source project support the study of the effect of surface roughness on the deformation behavior of flexible plates in flow fields？

Dear developers，
It is my honor to write this letter！ I hope to investigate the effect of surface roughness on the deformation behavior of flexible plates in flow fields. However, due to a lack of understanding of the underlying logic of the project, I am not sure whether this open project can achieve my expected goals and whether there is complete validation of reference data and literature. I hope to receive a response from the developers, and I will be very grateful! Thank you!

The result doesn't deform in the visualization

I run the demo problems: turtle_demo and TF_fsi. The result doesn't deform as illustrated in the attached GIF on the readme page when I input them into ParaView. I think it may because I didn't input those files in the correct way.

Can someone tell me how should I import the XDMF and H5 files to ParaView?

posssibility to extent the fluid solver to be turbulent

Dear everyone

Recently I have been searching for a way to implement turbulent flow simulation in turtleFSI.
There is a post in fenics forum:
https://fenicsproject.discourse.group/t/navier-stokes-simulation-getting-slower-with-time/11452
He seems to have successfully implemented the LES Smagorinsky model in FENICSX.
Do you think it is possible to extend the existing fluid model to LES in turtlefsi?

Best regards
ZHANG

Found no facets matching domain for boundary condition.(TF_csm)

Dear everyone

When testing TF_csm, I tried to change the fixed boundary conditions for the solid domain by modifying the def get_mesh_domain_and_boundaries

def get_mesh_domain_and_boundaries(c_x, c_y, f_H, f_L, R, **namespace):
    # Read mesh
    mesh = Mesh(path.join(path.dirname(path.abspath(__file__)),  "mesh", "TF_csm.xml.gz"))
    mesh = refine(mesh)

    # Mark boundaries
    Barwall =  AutoSubDomain(lambda x: ((x[0] - c_x)**2 + (x[1] - c_y)**2 < R**2 + DOLFIN_EPS * 1e5))
    #Downwall = AutoSubDomain(lambda x: near(x[1], c_y - f_H / 2))#0.2-0.02/2=0.19
    #Upwall =   AutoSubDomain(lambda x: near(x[1], c_y + f_H / 2))#0.2+0.02/2=0.21
    #Rightwall= AutoSubDomain(lambda x: near(x[0], c_x + R + f_L))#0.2+0.05+0.35=0.60
 
    def Down(x, on_boundary):
        return near(x[1], 0.19) and on_boundary
    def Up(x, on_boundary):
        return near(x[1], 0.21) and on_boundary
    def Right(x, on_boundary):
        return near(x[0], 0.60) and on_boundary

    Downwall = AutoSubDomain(Down)
    Upwall= AutoSubDomain(Up)
    Rightwall= AutoSubDomain(Right)
    
    boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)
    boundaries.set_all(0)
    Barwall.mark(boundaries, 1)
    Downwall.mark(boundaries, 2)
    Upwall.mark(boundaries, 3)
    Rightwall.mark(boundaries, 4)

    # Mark domain
    domains = MeshFunction("size_t", mesh, mesh.geometry().dim())
    domains.set_all(1)

    return mesh, domains, boundaries

def create_bcs(DVP, boundaries, **namespace):
    # Clamp on the right hand side
    u_barwall = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 4)

    return dict(bcs=[u_barwall])

But it seems it can't be found during calculating.

hwloc/linux: Ignoring PCI device with non-16bit domain.
Pass --enable-32bits-pci-domain to configure to support such devices
(warning: it would break the library ABI, don't enable unless really needed).
*** Warning: Found no facets matching domain for boundary condition.
*** Warning: Found no facets matching domain for boundary condition.
*** Warning: Found no facets matching domain for boundary condition.
*** Warning: Found no facets matching domain for boundary condition.
WARNING: Could not import module "cppimport" - it is required for using
 Probes, but not installed on your system. Install with "pip install cppimport"
This is used for solid using electic tensor C ----ZHANG
1 solid region(s) found, using following parameters
{'dx_s_id': 1, 'material_model': 'C_arbitrary', 'rho_s': 5550, 'c11': 226000000000.0, 'c13': 124000000000.0, 'c33': 216000000000.0, 'c44': 44200000000.0}
1 fluid region(s) found, using following parameters
{'dx_f_id': 0, 'rho_f': 1000.0, 'mu_f': 1.0}
dx_f_id_list: [0] 
dx_s_id_list: [1] 
Using fluid file: no_fluid 
Using solid file: solid_C_arbitraty_reversetest 
Using extrapolation file: no_extrapolation 
Compute Jacobian matrix
Newton iteration 0: r (atol) = 7.506e+00 (tol = 1.000e-07), r (rel) = 8.203e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 1.434e-05 (tol = 1.000e-07), r (rel) = 5.072e-09 (tol = 1.000e-07) 
Distance x: 2.852964e-18
Distance y: -4.900000e-04
Solved for timestep 1, t = 0.0100 in 0.4 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 7.506e+05 (tol = 1.000e-07), r (rel) = 8.204e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 4.308e-05 (tol = 1.000e-07), r (rel) = 1.521e-08 (tol = 1.000e-07) 
Distance x: 1.569242e-17
Distance y: -1.960000e-03
Solved for timestep 2, t = 0.0200 in 0.3 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 1.501e+06 (tol = 1.000e-07), r (rel) = 8.206e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 7.208e-05 (tol = 1.000e-07), r (rel) = 2.536e-08 (tol = 1.000e-07) 
Distance x: 4.674880e-17
Distance y: -4.410000e-03
Solved for timestep 3, t = 0.0300 in 0.3 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 2.252e+06 (tol = 1.000e-07), r (rel) = 8.208e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 1.017e-04 (tol = 1.000e-07), r (rel) = 3.550e-08 (tol = 1.000e-07) 
Distance x: 1.055812e-16
Distance y: -7.840000e-03
Solved for timestep 4, t = 0.0400 in 0.4 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 3.002e+06 (tol = 1.000e-07), r (rel) = 8.212e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 1.302e-04 (tol = 1.000e-07), r (rel) = 4.564e-08 (tol = 1.000e-07) 
Distance x: 1.982616e-16
Distance y: -1.225000e-02
Solved for timestep 5, t = 0.0500 in 0.3 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 3.753e+06 (tol = 1.000e-07), r (rel) = 8.216e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 1.614e-04 (tol = 1.000e-07), r (rel) = 5.579e-08 (tol = 1.000e-07) 
Distance x: 3.371549e-16
Distance y: -1.764000e-02
Solved for timestep 6, t = 0.0600 in 0.3 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 4.504e+06 (tol = 1.000e-07), r (rel) = 8.221e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 1.905e-04 (tol = 1.000e-07), r (rel) = 6.593e-08 (tol = 1.000e-07) 
Distance x: 5.248376e-16
Distance y: -2.401000e-02
Solved for timestep 7, t = 0.0700 in 0.3 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 5.254e+06 (tol = 1.000e-07), r (rel) = 8.226e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 2.260e-04 (tol = 1.000e-07), r (rel) = 7.607e-08 (tol = 1.000e-07) 
Distance x: 7.799181e-16
Distance y: -3.136000e-02
Solved for timestep 8, t = 0.0800 in 0.3 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 6.005e+06 (tol = 1.000e-07), r (rel) = 8.233e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 2.565e-04 (tol = 1.000e-07), r (rel) = 8.622e-08 (tol = 1.000e-07) 
Distance x: 1.094141e-15
Distance y: -3.969000e-02
Solved for timestep 9, t = 0.0900 in 0.3 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 6.756e+06 (tol = 1.000e-07), r (rel) = 8.240e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 2.872e-04 (tol = 1.000e-07), r (rel) = 9.636e-08 (tol = 1.000e-07) 
Distance x: 1.496842e-15
Distance y: -4.900000e-02
Solved for timestep 10, t = 0.1000 in 0.4 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 7.506e+06 (tol = 1.000e-07), r (rel) = 8.248e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 3.184e-04 (tol = 1.000e-07), r (rel) = 1.065e-07 (tol = 1.000e-07) 
Newton iteration 2: r (atol) = 1.042e-04 (tol = 1.000e-07), r (rel) = 9.129e-14 (tol = 1.000e-07) 
Distance x: 1.939238e-15
Distance y: -5.929000e-02
Solved for timestep 11, t = 0.1100 in 0.4 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 8.257e+06 (tol = 1.000e-07), r (rel) = 8.257e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 3.611e-04 (tol = 1.000e-07), r (rel) = 1.166e-07 (tol = 1.000e-07) 
Newton iteration 2: r (atol) = 1.484e-04 (tol = 1.000e-07), r (rel) = 9.999e-14 (tol = 1.000e-07) 
Distance x: 2.379777e-15
Distance y: -7.056000e-02
Solved for timestep 12, t = 0.1200 in 0.5 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 9.007e+06 (tol = 1.000e-07), r (rel) = 8.267e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 3.894e-04 (tol = 1.000e-07), r (rel) = 1.268e-07 (tol = 1.000e-07) 
Newton iteration 2: r (atol) = 1.745e-04 (tol = 1.000e-07), r (rel) = 1.087e-13 (tol = 1.000e-07) 
Distance x: 2.823974e-15
Distance y: -8.281000e-02
Solved for timestep 13, t = 0.1300 in 0.5 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 9.758e+06 (tol = 1.000e-07), r (rel) = 8.278e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 4.276e-04 (tol = 1.000e-07), r (rel) = 1.369e-07 (tol = 1.000e-07) 
Newton iteration 2: r (atol) = 1.915e-04 (tol = 1.000e-07), r (rel) = 1.174e-13 (tol = 1.000e-07) 
Distance x: 3.263889e-15
Distance y: -9.604000e-02
Solved for timestep 14, t = 0.1400 in 0.4 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 1.051e+07 (tol = 1.000e-07), r (rel) = 8.289e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 4.574e-04 (tol = 1.000e-07), r (rel) = 1.471e-07 (tol = 1.000e-07) 
Newton iteration 2: r (atol) = 1.915e-04 (tol = 1.000e-07), r (rel) = 1.261e-13 (tol = 1.000e-07) 
Distance x: 3.721201e-15
Distance y: -1.102500e-01
Solved for timestep 15, t = 0.1500 in 0.3 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 1.126e+07 (tol = 1.000e-07), r (rel) = 8.301e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 5.223e-04 (tol = 1.000e-07), r (rel) = 1.572e-07 (tol = 1.000e-07) 
Newton iteration 2: r (atol) = 2.812e-04 (tol = 1.000e-07), r (rel) = 1.348e-13 (tol = 1.000e-07) 
Distance x: 4.157062e-15
Distance y: -1.254400e-01
Solved for timestep 16, t = 0.1600 in 0.4 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 1.201e+07 (tol = 1.000e-07), r (rel) = 8.314e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 5.593e-04 (tol = 1.000e-07), r (rel) = 1.674e-07 (tol = 1.000e-07) 
Newton iteration 2: r (atol) = 2.982e-04 (tol = 1.000e-07), r (rel) = 1.435e-13 (tol = 1.000e-07) 
Distance x: 4.616574e-15
Distance y: -1.416100e-01
Solved for timestep 17, t = 0.1700 in 0.3 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 1.276e+07 (tol = 1.000e-07), r (rel) = 8.328e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 5.893e-04 (tol = 1.000e-07), r (rel) = 1.775e-07 (tol = 1.000e-07) 
Newton iteration 2: r (atol) = 3.199e-04 (tol = 1.000e-07), r (rel) = 1.522e-13 (tol = 1.000e-07) 
Distance x: 5.059804e-15
Distance y: -1.587600e-01
Solved for timestep 18, t = 0.1800 in 0.3 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 1.351e+07 (tol = 1.000e-07), r (rel) = 8.343e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 6.376e-04 (tol = 1.000e-07), r (rel) = 1.876e-07 (tol = 1.000e-07) 
Newton iteration 2: r (atol) = 3.783e-04 (tol = 1.000e-07), r (rel) = 1.608e-13 (tol = 1.000e-07) 
Distance x: 5.505911e-15
Distance y: -1.768900e-01
Solved for timestep 19, t = 0.1900 in 0.3 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 1.426e+07 (tol = 1.000e-07), r (rel) = 8.358e-03 (tol = 1.000e-07) 
Newton iteration 1: r (atol) = 6.766e-04 (tol = 1.000e-07), r (rel) = 1.978e-07 (tol = 1.000e-07) 
Newton iteration 2: r (atol) = 3.762e-04 (tol = 1.000e-07), r (rel) = 1.695e-13 (tol = 1.000e-07) 
Distance x: 5.981709e-15
Distance y: -1.960000e-01
Solved for timestep 20, t = 0.2000 in 0.6 s
Compute Jacobian matrix
Newton iteration 0: r (atol) = 1.501e+07 (tol = 1.000e-07), r (rel) = 8.374e-03 (tol = 1.000e-07) 
*** Warning: Found no facets matching domain for boundary condition.
*** Warning: Found no facets matching domain for boundary condition.
*** Warning: Found no facets matching domain for boundary condition.
*** Warning: Found no facets matching domain for boundary condition.
*** Warning: Found no facets matching domain for boundary condition.
*** Warning: Found no facets matching domain for boundary condition.

And actually, in the four boundaries, both top and right can't be found but left and bottom seem work fine. Is there anything wrong for the boundary conditions I defined?

Pressure should not be part of save_deg >1 ?

turtleFSI/turtleFSI/problems/__init__.py

Lines 235 to 247 in d149375

	pe = FiniteElement('CG', mesh.ufl_cell(), p_deg)
	FSdv = FunctionSpace(mesh, dve) # Higher degree FunctionSpace for d and v
	FSp= FunctionSpace(mesh, pe) # Higher degree FunctionSpace for p
	# Copy mesh
	mesh_viz = Mesh(mesh)
	for i in range(save_deg-1):
	mesh_viz = refine(mesh_viz) # refine the mesh
	# Create visualization function space for d, v and p
	dve_viz = VectorElement('CG', mesh_viz.ufl_cell(), 1)
	pe_viz = FiniteElement('CG', mesh_viz.ufl_cell(), 1)

Since we use p_deg=1 for most of the case, it does not make sense that the pressure is part of save_deg >1 functionality.
We should check if the original degree is actually higher than 1st order before calling save_deg>1 part. Otherwise, we are interpolating the P1 solution from coarse to fine mesh.

Fix 1st lame coefficient in `TF_fsi`

turtleFSI/turtleFSI/problems/TF_fsi.py

Line 41 in e1b5350

lambda_s=4e6, # Solid 1st Lame Coefficient [Pa]

should be lambda_s=8e6

merge_visualization_files not called on restart

merge_visualization_files is not called on restart because restart_folder is set to None in monolithic.py#L159.

checkpoint is only saving one time step, n-1

checkpoint function states that it saves last two time steps, but it only saves one time step.

turtleFSI/turtleFSI/problems/__init__.py

Lines 183 to 200 in e8574bb

	def checkpoint(dvp_, default_variables, checkpoint_folder, mesh, **namespace):
	"""Utility function for storing the current parameters and the last two time steps"""
	# Only update variables that exists in default_variables
	default_variables.update((k, namespace[k]) for k in (default_variables.keys() & namespace.keys()))
	# Dump default parameters
	if MPI.rank(MPI.comm_world) == 0:
	with open(str(checkpoint_folder.joinpath("default_variables.pickle")), "bw") as f:
	pickle.dump(default_variables, f)
	# Dump physical fields
	fields = _get_fields(dvp_, mesh)
	# Write fields to temporary file to avoid corruption of existing checkpoint
	for name, field in fields:
	checkpoint_path = str(checkpoint_folder.joinpath("tmp_" + name + ".xdmf"))
	with XDMFFile(MPI.comm_world, checkpoint_path) as f:
	f.write_checkpoint(field, name)

field is generated inside

turtleFSI/turtleFSI/problems/__init__.py

Lines 364 to 374 in e8574bb

	def _get_fields(dvp_, mesh):
	d1 = dvp_["n-1"].sub(0, deepcopy=True)
	v1 = dvp_["n-1"].sub(1, deepcopy=True)
	p1 = dvp_["n-1"].sub(2, deepcopy=True)
	fields = [('d1', d1), ('v1', v1), ('p1', p1)]
	if len(dvp_["n-1"]) == mesh.geometric_dimension() * 3 + 1:
	w1 = dvp_["n-1"].sub(3, deepcopy=True)
	fields += [('w1', w1)]
	return fields

and you can see that it only has dvp_["n-1”]. So when restarting the simulation, we are right now assigning n-1 to both n and n-1

turtleFSI/turtleFSI/problems/__init__.py

Lines 356 to 361 in e8574bb

	assign(dvp_["n-1"].sub(0), fields[0][1]) # update d_["n-1"] to checkpoint d_["n-1"]
	assign(dvp_["n-1"].sub(1), fields[1][1]) # update v_["n-1"] to checkpoint v_["n-1"]
	assign(dvp_["n-1"].sub(2), fields[2][1]) # update p_["n-1"] to checkpoint p_["n-1"]
	assign(dvp_["n"].sub(0), fields[0][1]) # update d_["n-1"] to checkpoint d_["n-1"]
	assign(dvp_["n"].sub(1), fields[1][1]) # update v_["n-1"] to checkpoint v_["n-1"]
	assign(dvp_["n"].sub(2), fields[2][1]) # update p_["n-1"] to checkpoint p_["n-1"]

This restart need to be fixed so that we actually restart the problem using two last time steps.

`configparser` convert all key names to lowercase

I noticed that configparser convert all key names to lowercase by its default. This is not a big problem but is not mentioned anywhere in our example and we actually use parameter with uppercase. As long as there is no demand for supporting upper case parameter, I will change example.config and turtle_demo to make it consistent.

Ref.

turtleFSI/docs/examples/example.config

Line 50 in e869c8d

Um=0.8

This will be treated as um

Replace `up_solve` with petsc KSP

Should be easier to access debug mode: https://github.com/KVSlab/turtleFSI/search?q=up_sol

right hand side vector `b` newly created at every time step?

I noticed that right hand side vector b is newly created in every time step as it is not returned at the end of newton iteration.
Here is the print of ID from first five time step

ID of b: 4335405648 at time step 0
ID of b: 6296605200 at time step 1
ID of b: 6309206976 at time step 2
ID of b: 6311653712 at time step 3
ID of b: 6305607536 at time step 4
ID of b: 6413148672 at time step 5

Since b is originally set to None inside solver_setup function, b becomes None again after completing one time step, and turtle assign new b in every time step as b is None at the beginning of every time step. Specifically, line 85 of newtonsolver.py gets called at the beginning of every time step instead of line 87

turtleFSI/turtleFSI/modules/newtonsolver.py

Lines 84 to 87 in 34925be

	if b is None:
	b = assemble(-F)
	else:
	assemble(-F, tensor=b)

I do not know how much impact it has, but I assume it is a waste of memory to create new b. Just adding b as a part of return arguments at the end of newton iteration should fix the problem.

turtleFSI/turtleFSI/modules/newtonsolver.py

Line 114 in 34925be

return dict(up_sol=up_sol, A=A)

- return dict(up_sol=up_sol, A=A)
+ return dict(up_sol=up_sol, A=A, b=b)

`merge_xml_files` does not work above python3.9

Since getchildren() is deprecated from python3.9, merge_xml_files does not work above python3.9

turtleFSI/turtleFSI/problems/__init__.py

Line 422 in b4df373

ind = 1 if len(last_node.getchildren()) == 3 else 2

I modified the function and it works with python3.10 on my laptop, but I do not know if it works with python3.8 or lower.
Can @johannesring or @jorgensd maybe test if replacing the merge_xml_files as follows works fine with python3.8?
You can use Taylor green problem that I mentioned in the other issue to test merge_xml_files.
Here is the link to the branch https://github.com/keiyamamo/turtleFSI/tree/fix_save_viz
and also new mega_xml_file (in case you just want to copy this function).

def merge_xml_files(files):
    # Get first timestep and trees
    first_timesteps = []
    trees = []
    for f in files:
        trees.append(ET.parse(f))
        root = trees[-1].getroot()
        first_timesteps.append(float(root[0][0][0][2].attrib["Value"]))

    # Index valued sort (bypass numpy dependency)
    first_timestep_sorted = sorted(first_timesteps)
    indexes = [first_timesteps.index(i) for i in first_timestep_sorted]

    # Get last timestep of first tree
    base_tree = trees[indexes[0]]
    last_node = base_tree.getroot()[0][0][-1]
    ind = 2 if len(last_node.findall('*')) == 4 else 1
    last_timestep = float(last_node[ind].attrib["Value"])

    # Append
    for index in indexes[1:]:
        tree = trees[index]
        for node in tree.getroot()[0][0].findall('*'):
            ind = 2 if len(node.findall('*')) == 4 else 1
            if last_timestep < float(node[ind].attrib["Value"]):
                base_tree.getroot()[0][0].append(node)
                last_timestep = float(node[ind].attrib["Value"])

    # Seperate xdmf files
    new_file = [f for f in files if "_" not in f.name.__str__()]
    old_files = [f for f in files if "_" in f.name.__str__()]

    # Write new xdmf file
    base_tree.write(new_file[0], xml_declaration=True)

    # Delete xdmf file
    [f.unlink() for f in old_files]

Adding multi-solid regions

Dear everyone

I meet a problem when adding the multi-solid regions with different solid materials by modifying the fsi3 case.

The tag of boundaries and regions set in gmsh are:

1 10 "Solidinterface"
1 11 "Inlet"
1 12 "Outlet"
1 13 "Wall"
1 14 "Bar"
1 15 "Barwall"
1 16 "Circle"
2 17 "Fluid" #region
2 18 "Solid1" #region
2 19 "Solid2" #region

To do it, I referred to the turtleFSI-master/turtleFSI/modules/domains.py and solid.py and changed them into domain_solid.py and solid_multiphase.py as follows:

#domain_solid.py
from turtleFSI.modules import *
from dolfin import Constant, inner, inv, grad, div

def assign_solid_domain_properties(dx, dx_s_id, rho_s, nu_s, mu_s, lambda_s, solid_properties, domains, **namespace):
    """
    Assigns solid properties to each region.   
    """

    # 1. Create differential for each solid region, and organize into solid_properties list of dicts
    dx_s = {}
    if isinstance(dx_s_id, list): # If dx_s_id is a list (i.e, if there are multiple solid regions):
        for solid_region in range(len(dx_s_id)):
            dx_s[solid_region] = dx(dx_s_id[solid_region], subdomain_data=domains) # Create dx_s for each solid region
        dx_s_id_list=dx_s_id
    else:
        dx_s[0] = dx(dx_s_id, subdomain_data=domains)
        dx_s_id_list=[dx_s_id]



    # Assign material properties to each region
    if len(solid_properties) == 0:
        if isinstance(dx_s_id, list): # If dx_s_id is a list (i.e, if there are multiple solid regions):
            for solid_region in range(len(dx_s_id)):
                solid_properties.append(
                    {"dx_s_id": dx_s_id[solid_region],
                     "rho_s": rho_s[solid_region],
                     "nu_s": nu_s[solid_region],
                     "mu_s": mu_s[solid_region],
                     "lambda_s": lambda_s[solid_region]})
        else:
            solid_properties.append(
                {"dx_s_id": dx_s_id,
                 "rho_s": rho_s,
                 "nu_s": nu_s,
                 "mu_s": mu_s,
                 "lambda_s": lambda_s})
    elif isinstance(solid_properties, dict):
        solid_properties = [solid_properties]


    return dict(dx_s=dx_s, dx_s_id_list=dx_s_id_list, solid_properties=solid_properties)

and

#solid_multiphase.py
from turtleFSI.modules import *
from dolfin import Constant, inner, grad

def solid_setup(d_, v_, phi, psi, dx_s, dx_s_id_list, solid_properties, k, theta,
                gravity, **namespace):

    delta = 1E7

    # Theta scheme constants
    theta0 = Constant(theta)
    theta1 = Constant(1 - theta)


    F_solid_linear = 0
    F_solid_nonlinear = 0
    for solid_region in range(len(dx_s_id_list)):
        rho_s = solid_properties[solid_region]["rho_s"]
        nu_s  = solid_properties[solid_region]["nu_s"]
        mu_s  = solid_properties[solid_region]["mu_s"]
        lambda_s = solid_properties[solid_region]["lambda_s"]

        # Temporal term and convection
        F_solid_linear = (rho_s/k * inner(v_["n"] - v_["n-1"], psi)*dx_s[solid_region]
                          + delta * rho_s * (1 / k) * inner(d_["n"] - d_["n-1"], phi) * dx_s[solid_region]
                          - delta * rho_s * inner(theta0 * v_["n"] + theta1 * v_["n-1"], phi) * dx_s[solid_region])

        # Gravity
        if gravity is not None:
            F_solid_linear -= inner(Constant((0, -gravity * rho_s)), psi)*dx_s[solid_region]

        # Stress
        F_solid_nonlinear = theta0 * inner(Piola1(d_["n"], lambda_s, mu_s), grad(psi)) * dx_s[solid_region]
        F_solid_linear += theta1 * inner(Piola1(d_["n-1"], lambda_s, mu_s), grad(psi)) * dx_s[solid_region]

    return dict(F_solid_linear=F_solid_linear, F_solid_nonlinear=F_solid_nonlinear)

The input parameters are

        # Temporal variables
        T=10,                         # End time [s]
        dt=0.001,                      # Time step [s]
        theta=0.51,                    # Temporal scheme

        # Physical constants ('FSI 3')
        Um=2.0,                       # Max. velocity inlet, CDF3: 2.0 [m/s]
        rho_f=1.0e3,                  # Fluid density [kg/m3]
        mu_f=1.0,                     # Fluid dynamic viscosity [Pa.s]

        solid_properties=[{"dx_s_id": 18, "rho_s": 2.0E3, "nu_s": 0.8, "mu_s": 4.0e6, "lambda_s": 8E6}, {"dx_s_id": 19, "rho_s": 1.0E3, "nu_s": 0.4, "mu_s": 2.0e6, "lambda_s": 4e6}],

        # Domain settings
        dx_f_id=17,       # Domain id of the fluid domain
        dx_s_id=[18,19],  # Domain id of the solid domain
        fluid="fluid",                             # ["fluid", "no_fluid"] Turn off fluid and only solve the solid problem
        solid="solid_multiphase",                  # ["solid", "no_solid"] Turn off solid and only solve the fluid problem

        # Problem specific
        folder="CylinderFlap_multisolidphase_interface_results",      # Name of the results folder
        extrapolation="biharmonic",   # No displacement to extrapolate
        extrapolation_sub_type="constrained_disp",  # Biharmonic type

But the results are bad:

It seems the code cannot rightly recognize the physical region, do some friends here know why? Is that because of the interface between two solid regions?

Thanks in advance!

Best regards
ZHANG

Update pytest

Current pytest implementation does not test mpirun or restart function. These two must be added to the pytest.
Also, I think we should reconsider TF_fsi as a benchmark problem and move to different problem with analytical solution.

xdmf file not correct when using `save_deg=2` and restart with different mesh patition

I found that xdmf files become ‘ill’ when we use save_deg=2 and also restart with different mesh partition.
I create a problem file that everyone can easily test and see how it looks like when this problem occurs.

Here is the link to the 2D Taylor-Green problem that’s designed to stop after 10 seconds.
https://github.com/keiyamamo/turtleFSI/blob/0208e536e4edce3efffcc450865896fe95a25633/turtleFSI/problems/tg2d.py

First, you run the problem with just single processor

 python monolithic.py -p tg2d

and then turtle will automatically stop after 10 seconds with the following message (the actual message may differ depending on your machine)

Solved for timestep 17, t = 0.1700 in 0.6 s
Given killtime reached! Stopping simulations cleanly...
Total simulation time 10.046949

After that you restart the problem with 2 processors with

mpirun -np 2 python monolithic.py -p tg2d --restart-folder tg2d_results/1/

this gives

Solved for timestep 18, t = 0.1800 in 0.5 s
:
:
Solved for timestep 46, t = 0.4600 in 0.4 s
Given killtime reached! Stopping simulations cleanly...
Total simulation time 10.176421
Merging visualization files

meaning that restart started from timestep 18 and ended at timestep 46 and merged visualization files.
However, if you open the xdmf file you see that the result does not look correct after timestep 18

This is because xdmf file points to the mesh information from velocity_run_1.h5 while velocity.h5 and velocity_run_1.h5 have different mesh information.

If you open, velocity.xdmf, you can see it points to velocity_run_1.h5:/Mesh/0/mesh/geometry

<Grid Name="mesh" GridType="Uniform">
        <Topology NumberOfElements="3200" TopologyType="Triangle" NodesPerElement="3">
          <DataItem Dimensions="3200 3" NumberType="UInt" Format="HDF">velocity_run_1.h5:/Mesh/0/mesh/topology</DataItem>
        </Topology>
        <Geometry GeometryType="XY">
          <DataItem Dimensions="1681 2" Format="HDF">velocity_run_1.h5:/Mesh/0/mesh/geometry</DataItem>
        </Geometry>

accessing external files on docker

Im currently using turtle fsi on docker and microsoft visual studio and have trouble accessing external filed and when i execute moving cylinder i get fenics@df3ab8a8010e:~/local/turtleFSI$ turtleFSI --problem MovingCylinder
{'A_ratio': 0.25,
'D': 0.1,
'F_r': 1.0,
'Re': 500,
'St': 0.228,
'T': 5,
'Um': 0.8,
'atol': 1e-07,
'bc_ids': [],
'c_s': 0.0,
'checkpoint_step': 500,
'compiler_parameters': {'add_tabulate_tensor_timing': False,
'cache_dir': '',
'convert_exceptions_to_warnings': False,
'cpp_optimize': True,
'cpp_optimize_flags': '-O2',
'epsilon': 1e-14,
'error_control': False,
'external_include_dirs': '',
'external_includes': '',
'external_libraries': '',
'external_library_dirs': '',
'form_postfix': False,
'format': 'ufc',
'generate_dummy_tabulate_tensor': False,
'log_level': 25,
'log_prefix': '',
'max_signature_length': 0,
'no-evaluate_basis_derivatives': True,
'optimize': True,
'output_dir': '.',
'precision': None,
'quadrature_degree': 4,
'quadrature_rule': None,
'representation': 'auto',
'split': False},
'counter': 0,
'd_deg': 1,
'ds_s_id': None,
'dt': 0.000125,
'dx_f_id': 1,
'dx_s_id': 2,
'extrapolation': 'laplace',
'extrapolation_sub_type': 'constant',
'factor': 50.0,
'fluid': 'fluid',
'fluid_properties': [],
'folder': 'results_moving_cylinder',
'gravity': None,
'k_s': 0.0,
'killtime': None,
'lambda_s': 450000.0,
'linear_solver': 'mumps',
'lmbda': 1.0,
'loglevel': 20,
'material_model': 'StVenantKirchoff',
'max_it': 50,
'mu_f': 0.2,
'mu_s': 50000.0,
'nu_s': 0.45,
'p_deg': 1,
'problem': 'MovingCylinder',
'recompute': 25,
'recompute_tstep': 100,
'restart_folder': None,
'rho_f': 1000,
'rho_s': 1000.0,
'robin_bc': False,
'rtol': 1e-07,
'save_deg': 1,
'save_step': 100,
'solid': 'no_solid',
'solid_properties': [],
'solver': 'newtonsolver',
'sub_folder': None,
't': 0,
'theta': 0.500125,
'u_inf': 1.0,
'v_deg': 1,
'verbose': True}
Traceback (most recent call last):
File "/usr/local/bin/turtleFSI", line 11, in
load_entry_point('turtleFSI', 'console_scripts', 'turtleFSI')()
File "/home/fenics/local/turtleFSI/turtleFSI/run_turtle.py", line 18, in main
from turtleFSI import monolithic
File "/home/fenics/local/turtleFSI/turtleFSI/monolithic.py", line 63, in
mesh, domains, boundaries = get_mesh_domain_and_boundaries(**vars())
File "/home/fenics/local/turtleFSI/turtleFSI/problems/MovingCylinder.py", line 81, in get_mesh_domain_and_boundaries
infile.read(mesh)
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

---

*** [email protected]

---

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

---

*** -------------------------------------------------------------------------
*** Error: Unable to open XDMF file.
*** Reason: XDMF file "mesh/MovingCylinder/mesh4.xdmf" does not exist.
*** Where: This error was encountered inside XDMFFile.cpp.
*** Process: 0

---

*** DOLFIN version: 2019.1.0
*** Git changeset: 74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

fenics@df3ab8a8010e:~/local/turtleFSI$

Help me please

How to get the same post results as the example case shown in turtleFSI

Thanks for the great work on this solver.

I ran the example cases of TF_fsi, and it ran and finished without any errors, but seeing the xdmf files in the 'TF_fsi_results\1\Visualization' I think the velocity contour is not the fsi result, because I didn't see any mesh deformation in the contour.

I think the demo fig is right, and it is really different from my contour. This is really strange.

But when I check the output x and y displacements and lift and drag, they are well presented.

May I know how to get the right contour shown in this website?
Which software do you use?
Is my test simulation wrong or is it just the way of postprocessing?

Very slow Newton iteration

I’m experiencing very slow Newton iteration due to slight increase in the absolute tolerance during the Newton iteration. In addition, there is a repetitive behavior of the atol and rel, meaning there are basically two or three values repeating after recomputing Jacobian. Following is the example from my simulation

Newton iteration 0: r (atol) = 1.587e-03 (tol = 1.000e-06), r (rel) = 2.483e-03 (tol = 1.000e-06) 
Newton iteration 1: r (atol) = 2.033e-05 (tol = 1.000e-06), r (rel) = 2.642e-03 (tol = 1.000e-06) 
Newton iteration 2: r (atol) = 2.043e-05 (tol = 1.000e-06), r (rel) = 4.638e-05 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 3: r (atol) = 2.043e-05 (tol = 1.000e-06), r (rel) = 1.067e-04 (tol = 1.000e-06) 
Newton iteration 4: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 8.132e-06 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 5: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 9.712e-05 (tol = 1.000e-06) 
Newton iteration 6: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 8.698e-06 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 7: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 9.601e-05 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 8: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.036e-04 (tol = 1.000e-06) 
Newton iteration 9: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 8.316e-06 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 10: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 9.582e-05 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 11: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.036e-04 (tol = 1.000e-06) 
Newton iteration 12: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 8.334e-06 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 13: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 9.582e-05 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 14: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.036e-04 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 15: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 1.032e-04 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 16: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.031e-04 (tol = 1.000e-06) 
Newton iteration 17: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 8.275e-06 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 18: r (atol) = 2.043e-05 (tol = 1.000e-06), r (rel) = 9.581e-05 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 19: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.036e-04 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 20: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 1.032e-04 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 21: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.031e-04 (tol = 1.000e-06) 
Newton iteration 22: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 8.275e-06 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 23: r (atol) = 2.043e-05 (tol = 1.000e-06), r (rel) = 9.581e-05 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 24: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.036e-04 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 25: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 1.032e-04 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 26: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.031e-04 (tol = 1.000e-06) 
Newton iteration 27: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 8.275e-06 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 28: r (atol) = 2.043e-05 (tol = 1.000e-06), r (rel) = 9.581e-05 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 29: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.036e-04 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 30: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 1.032e-04 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 31: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.031e-04 (tol = 1.000e-06) 
Newton iteration 32: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 8.275e-06 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 33: r (atol) = 2.043e-05 (tol = 1.000e-06), r (rel) = 9.581e-05 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 34: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.036e-04 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 35: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 1.032e-04 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 36: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.031e-04 (tol = 1.000e-06) 
Newton iteration 37: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 8.275e-06 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 38: r (atol) = 2.043e-05 (tol = 1.000e-06), r (rel) = 9.581e-05 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 39: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.036e-04 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 40: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 1.032e-04 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 41: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.031e-04 (tol = 1.000e-06) 
Newton iteration 42: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 8.275e-06 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 43: r (atol) = 2.043e-05 (tol = 1.000e-06), r (rel) = 9.581e-05 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 44: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.036e-04 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 45: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 1.032e-04 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 46: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.031e-04 (tol = 1.000e-06) 
Newton iteration 47: r (atol) = 2.042e-05 (tol = 1.000e-06), r (rel) = 8.275e-06 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 48: r (atol) = 2.043e-05 (tol = 1.000e-06), r (rel) = 9.581e-05 (tol = 1.000e-06) 
Compute Jacobian matrix
Newton iteration 49: r (atol) = 2.040e-05 (tol = 1.000e-06), r (rel) = 1.036e-04 (tol = 1.000e-06) 
Solved for timestep 4950, t = 1.6823 in 2181.7 s

Notice that there are basically three values of atol and rel after the second iteration, and due to slight increase in the absolute tolerance (not the relative tolerance because we removed the test), turtle computes Jacobian almost every iteration. This results in 2181.7s, more than 35min, for just one time step. I would not say this is a bug in turtle, but might be something we could avoid? For example, does it converge better if we allow slight increase in the absolute tolerance instead of recomputing Jacobian? This problem might be just unlucky, but I think it’s worth sharing here. Could @jorgensd perhaps comment on this ? I guess I’m hitting some sort of local minimum but I’m not knowledgeable enough to say if this is healthy or unhealthy behaviour of turtle.

A question about CSM case

Dear everyone

When studying the Demo case TF_CSM, I got a question about the result.
As this case simulated an elastic beam adding only the gravitational force, I think the beam will start to deform and vibrate initially. However, the vibration amplitude will decrease and keep still at last.
Why the simulation results show it keeps vibrating?
I am not questioning the accuracy of this result because the reference paper also has this phenomenon. But I just don't know why the vibration will not decrease.

Best regards

write explicitly that `warp by vector` is necessary to see the mesh deformation

Because there were questions about how to see the mesh deformation, it would be nice to include it in the documentation

some functions can not be found in this version of fenics in turtlefsi

Dear everyone:

I want to develop my work on multi-physics problems and use the Lagrange multiplier to add some constraints to the boundaries.
So I tested a Poisson problem, and I attached its code below. However, when I use the turtlefsi provided fenics2019.1, it will show some error messages like:

    gamma_mesh = MeshView.create(marker, 1)
                 ^^^^^^^^
NameError: name 'MeshView' is not defined

I also noticed many other functions like subspace or mixedfunctionspace that cannot be used in this version of fenics2019.1. But when I reached the forum of fenics, they said it should be in the fenics 2019.1. That makes me curious.

Does someone know how to fix this problem? If I can use higher fenics and keep on using the turtlefsi? Or am I doing wrong with this provided fenics in turtlefsi?

from __future__ import print_function
from fenics import *
from dolfin import *
from ufl import nabla_div
import numpy

# Scaled variables
n = 10
mesh = UnitCubeMesh(n, n, n)
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
for f in facets(mesh):
    marker[f] = 0.5 - DOLFIN_EPS < f.midpoint().x() < 0.5 + DOLFIN_EPS
gamma_mesh = MeshView.create(marker, 1)

V = FunctionSpace(mesh, "CG", 1)
LM = FunctionSpace(gamma_mesh, "DG", 0)
W = MixedFunctionSpace(V, LM)
(u, l) = TrialFunctions(W)
(v, e) = TestFunctions(W)


def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS


bc = DirichletBC(V, Constant(0.0), boundary)

dV = Measure(" dx ", domain=mesh)
dL = Measure(" dx ", domain=gamma_mesh)
fu = Expression(" 10* exp ( -( pow (x [0] - 0.5 , 2) + pow ( x [1] - 0.5 , 2)) / 0.02) ")
fl = Constant(c)  # u = c on Gamma
a = inner(grad(u), grad(v)) * dV + v * l * dL + u * e * dL
L = fu * v * dV + fl * e * dL

# Solve
sol = Function(W)
solve(a == L, sol, bc, solver_parameters={"linear_solver": "direct"})

I appreciate any help you can provide.
Zhang

Replace `gravity` with `body force`

gravity should be replaced with body force. Current implementation only allows us to apply force in y-direction

turtleFSI/turtleFSI/modules/solid.py

Lines 63 to 66 in c52ebc2

	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]

Oasis can be used as a reference point for the body force.
https://github.com/mikaem/Oasis/blob/da712b87a8998d0669df0448b1a4bbd9fea284d6/oasis/NSfracStep.py#L158-L159

xdmf issues

Hi- I've been having trouble opening the xdmf files the code generates from the simulations in Paraview. This is in contrast to the mesh files provided, which load up just fine. I suppose the h5 files are problematic in some way? Looking around the web, it seems to be a common issue for a lot of people. Would you be able to suggest a workaround or a suitable alternative format to use? Thank you.

change from .pickle to .json?

default_variables.pickle should be in json format