Some questions in the solid.py

[zhangmuElias]

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

[keiyamamo]

Hi @zhangmuElias

For the choice of test function, we follow the formulation by this PhD thesis in which the structure problem is formulated as

In our case, we do not include damping term so $\gamma_{\omega}$ = $\gamma_s = 0$. Unfortunately, I cannot find the explanation about the choice of particular test function for each equation in the thesis, but I assume it is related to the stability or continuity necessary for the well-posedness of the problem. I think the choice of test functions come from functional analysis and not necessary easy to relate to the physics of the problem.

Because of the mixed formulation between elastic and continuity equation, I do not think you can have the same test function for both of the equations, and that’s why we have psi for the elastic equation while we use phi for the continuity equation. In our case, it might still work even if you switch phi and psi in all the formulations because we use the same finite element for the displacement and velocity, but I’m not sure if that’s worth investigating. I’m not sure if I answered you question, but hope it helps. @jorgensd can maybe give more correct explanation on this.

[zhangmuElias]

Dear @keiyamamo

Thanks again for answering me every time for my question. It really helped me a lot! There is one thing I want to make sure about what you said:

it might still work even if you switch phi and psi in all the formulations

Do you mean if the phi is used for the elastic equation and psi for the continuity equation, it will still go well for this code? What only matters is to use different test functions for the elastic and continuity equations?

Also, I can test this thought on my computer.

Best regards Runze ZHANG

[keiyamamo]

Hi @zhangmuElias

What is the motivation for you to switch test functions? My apology for the confusing my last reply, I honestly do not know if it will work or not. More accurately, I was thinking " (Switching test functions should not work if we had different finite element or polynomial order between displacement and velocity, but we use the same element with the same order on the same mesh, so ) it might still work."

[zhangmuElias]

Dear @keiyamamo

For the motivation, it's about the piezoelectric field I want to add to the solid part in turtlefsi. I will tell you in detail in next post. I just tried to switch the test functions:

    for solid_region in range(len(dx_s_id_list)):
        rho_s = solid_properties[solid_region]["rho_s"]

        ## Temporal term and convection 
        F_solid_linear += (rho_s/k * inner(v_["n"] - v_["n-1"], phi)*dx_s[solid_region]
                          + delta * rho_s * (1 / k) * inner(d_["n"] - d_["n-1"], psi) * dx_s[solid_region]
                          - delta * rho_s * inner(theta0 * v_["n"] + theta1 * v_["n-1"], psi) * 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_C_Arbitrary(d_["n"], solid_properties[solid_region]), (Identity(2) + grad(phi))) * dx_s[solid_region]
        F_solid_linear += theta1 * inner(Piola1_C_Arbitrary(d_["n-1"], solid_properties[solid_region]), (Identity(2) + grad(phi))) * dx_s[solid_region]

        # Gravity
        if gravity is not None and mesh.geometry().dim() == 2:
            F_solid_linear -= inner(Constant((0, -gravity * rho_s)), phi)*dx_s[solid_region]
        elif gravity is not None and mesh.geometry().dim() == 3:
            F_solid_linear -= inner(Constant((0, -gravity * rho_s,0)), phi)*dx_s[solid_region]

Because I am using csm case to test, so I only changed the solid.py file. And the errors came when I ran the solver:

Using fluid file: no_fluid 
Using solid file: solid_C_arbitraty_reversetest 
Using extrapolation file: no_extrapolation 
Calling FFC just-in-time (JIT) compiler, this may take some time.
Compute Jacobian matrix
Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/bin/turtleFSI_master_c", line 10, in <module>
    sys.exit(main())
             ^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/turtleFSI/run_turtle_master_c.py", line 18, in main
    from turtleFSI import monolithic_master_c
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/turtleFSI/monolithic_master_c.py", line 174, in <module>
    vars().update(newtonsolver(**vars()))
                  ^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/turtleFSI/modules_master/newtonsolver.py", line 67, in newtonsolver
    A = assemble(J_nonlinear, tensor=A,
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/dolfin/fem/assembling.py", line 198, in assemble
    dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/dolfin/fem/assembling.py", line 56, in _create_dolfin_form
    return Form(form,
           ^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/dolfin/fem/form.py", line 43, in __init__
    ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/dolfin/jit/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/dolfin/jit/jit.py", line 97, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/ffc/jitcompiler.py", line 130, in jit_build
    module, signature = dijitso.jit(jitable=ufl_object,
                        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/dijitso/jit.py", line 165, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
                                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/ffc/jitcompiler.py", line 65, in jit_generate
    code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
                                            ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/ffc/compiler.py", line 142, in compile_form
    return compile_ufl_objects(forms, "form", object_names,
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/ffc/compiler.py", line 185, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/ffc/analysis.py", line 89, in analyze_ufl_objects
    form_datas = tuple(_analyze_form(form, parameters)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/ffc/analysis.py", line 89, in <genexpr>
    form_datas = tuple(_analyze_form(form, parameters)
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/ffc/analysis.py", line 169, in _analyze_form
    form_data = compute_form_data(form,
                ^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/ufl/algorithms/compute_form_data.py", line 418, in compute_form_data
    check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)  # Currently testing how fast this is
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/ufl/algorithms/check_arities.py", line 177, in check_form_arity
    check_integrand_arity(itg.integrand(), arguments, complex_mode)
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/ufl/algorithms/check_arities.py", line 159, in check_integrand_arity
    arg_tuples = map_expr_dag(rules, expr, compress=False)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/ufl/corealg/map_dag.py", line 37, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/ufl/corealg/map_dag.py", line 86, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/ufl/algorithms/check_arities.py", line 48, in sum
    raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b)))
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments () vs ('v_0',).

Do you know why this happens?

[keiyamamo]

Have you changed no_fluid.py as well ? There is also psi and phi there although I'm not sure if that matters.

[zhangmuElias]

Yes, the original no_fluid.py writes:

        F_fluid_linear +=    inner(Constant(tuple([0] * mesh.geometry().dim())), psi) * dx_f[fluid_region]
        F_fluid_nonlinear += inner(Constant(tuple([0] * mesh.geometry().dim())), phi) * dx_f[fluid_region]

I changed it to:

        F_fluid_linear += inner(Constant(tuple([0] * mesh.geometry().dim())), phi) * dx_f[fluid_region]
        F_fluid_nonlinear += inner(Constant(tuple([0] * mesh.geometry().dim())), psi) * dx_f[fluid_region]

Actually, I tried both these two type, and they show the same errors.

[keiyamamo]

Hi @zhangmuElias

Please have a look at here or here. In turtleFSI, we do not really expect variational form to be changed, so unfortunately I cannot spend too much time on that.

[zhangmuElias]

As for what I am doing.

I already tried successfully coded the geometry linear (small deformation) solid problem with the piezoelectric field in fenics. Because I want to study the piezoelectric energy harvester in fsi, I am now trying to expand the piezoelectric field for turtlefsi.

I changed the solid.py into

    for solid_region in range(len(dx_s_id_list)):
        rho_s = solid_properties[solid_region]["rho_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]) # 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_piezo(d_["n"],φ_["n"], solid_properties[solid_region]), grad(psi)) * dx_s[solid_region]
        F_solid_linear += theta1 * inner(Piola1_piezo(d_["n-1"],φ_["n-1"], solid_properties[solid_region]), grad(psi)) * dx_s[solid_region]

        # piezo part non
        F_solid_linear +=  (delta *  (1 / k) * dot(φ_["n"] - φ_["n-1"], phi_e) * dx_s[solid_region]
                          - delta *  dot(theta0 * Φ_["n"] + theta1 * Φ_["n-1"], phi_e) * dx_s[solid_region])  # Maybe we could add viscoelasticity with v["n"] term instead of (1 / k) * inner(d_["n"] - d_["n-1"], phi)
        F_solid_nonlinear -= theta0* dot(e_disp(d_["n"],φ_["n"], solid_properties[solid_region]), -grad(psi_e)) * dx_s[solid_region]
        F_solid_linear -= theta1* dot(e_disp(d_["n-1"],φ_["n-1"], solid_properties[solid_region]), -grad(psi_e)) * 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] 
        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] 

Although when I ran my code, it didn't give errors but it diverged in the first time step. Compute Jacobian matrix

Newton iteration 21: r (atol) = 5.307e+07 (tol = 1.000e-07), r (rel) = 2.888e+18 (tol = 1.000e-07) 
Compute Jacobian matrix
Newton iteration 22: r (atol) = 6.191e+07 (tol = 1.000e-07), r (rel) = 1.283e+18 (tol = 1.000e-07) 
Compute Jacobian matrix
Newton iteration 23: r (atol) = 5.190e+07 (tol = 1.000e-07), r (rel) = 5.865e+18 (tol = 1.000e-07) 
Compute Jacobian matrix
Newton iteration 24: r (atol) = 9.419e+07 (tol = 1.000e-07), r (rel) = 2.264e+19 (tol = 1.000e-07) 
Compute Jacobian matrix
Newton iteration 25: r (atol) = 2.537e+08 (tol = 1.000e-07), r (rel) = 9.836e+17 (tol = 1.000e-07) 
Compute Jacobian matrix
Traceback (most recent call last):
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/bin/turtleFSI_master_piezo", line 10, in <module>
    sys.exit(main())
             ^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/turtleFSI/run_turtle_master_piezo.py", line 18, in main
    from turtleFSI import monolithic_master_piezo
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/turtleFSI/monolithic_master_piezo.py", line 180, in <module>
    vars().update(newtonsolver(**vars()))
                  ^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dell-123/anaconda3/envs/turtlefsiproject/lib/python3.11/site-packages/turtleFSI/modules_master/newtonsolver.py", line 92, in newtonsolver
    raise RuntimeError("Error: The simulation has diverged during the Newton solve.")
RuntimeError: Error: The simulation has diverged during the Newton solve. 

[keiyamamo]

This delta term could be specific to the problem and should not really be hard-coded in my opinion. If you are adding term to the continuity equation, you might have to tune the value, assuming your variational formulation is correct (which I cannot tell since I'm not familiar with your topic).

[zhangmuElias]

Dear @keiyamamo You have already helped a lot. Thank you! Ok, let me try different delta to see what will happen.

[keiyamamo]

Closing. Feel free to reopen if necessary @zhangmuElias