Skip to content

Comments

Add potential-based contact model#7

Open
LeiShi3374 wants to merge 79 commits intovvedula22:masterfrom
LeiShi3374:master
Open

Add potential-based contact model#7
LeiShi3374 wants to merge 79 commits intovvedula22:masterfrom
LeiShi3374:master

Conversation

@LeiShi3374
Copy link

Add potential-based contact model and modified the penalty-based contact model

aabrown100-git and others added 30 commits September 23, 2021 11:29
Created a new function based on GNNB called GNNBT (T for current/deformed configuration). This function returns a normal vector weighted by the Jacobian, where the normal vector is the surface normal in the current configuration, and the Jacobian is the Jacobian of the mapping from the parent surface element to the current configuration surface element. Using GNNBT instead of GNNB in IntegV means integrating over the current configuration surface, instead of over the reference configuration surface. In particular, this fixes the velocity flux calculation when using struct (I think ustruct too). Before, the velocity flux was calculated as an integral over the reference configuration surface, which is not accurate (at least when using the velocity flux as the rate of change of enclosed volume)
I went through the code to understand how to fix the cplBC tangent matrix contribution when follower pressure load is used. I think I need to modify how valM (in ADDBCMUL) is computed to take into account deformation
Adding more comments during my search for how cplBC adds to the tangent matrix. I believe the integrals in Moghadam et al. eq. 27 are computed in BAFINI.f in the variable sV. However, these are computed just once at initialization, and they are computed using the weighted normal vector in the reference configuration (using GNNB())
I added a piece of code to BNEUFOLWP() (for follower pressure load) taken from FSILSINI() which updates the variable lhs%val to take into account the deformed geometry. I believe this val variable is the integral of the shape function times the normal vector over a surface, which is exactly what is needed to compute the cplBC tangent contribution (Moghadam et al. eq. 27). I added this piece of code so that these integrals are recomputed every Newton iteration.
Just some variables undefined. Should be able to get it working soon.
Fixed the calculation of the resistance or coupled BC contribution to the tangent matrix when follower pressure load is used. When follower pressure load is used, the resistance BC contribution to the residual vector involves an integral over the deformed configuration surface. The contribution to the stiffness matrix used the integral over the reference configuration surface, even with follower pressure load on, likely because the resistance/coupled BC code was developed for fixed domain CFD. As a result, the resistance BC contribution to the stiffness matrix was inconsistent with the contribution to the residual. Before this change, the struct-LPN had poor convergence when the deformed config was significantly different from the reference config. With this change, convergence is excellent even with large deformation.
Adding comments explaining in greater detail the changes I have made recently to svFSI, relating specifically to

1) changing the volume flux calculation for struct to account for deformed geometry

2) changing the calculation of the integrals involved in the resistance BC contribution to the tangent matrix.
Made changes to NN.f so that a cap mesh, which has surface element that do not lie on volume elements, is still processed and an area weighted normal vector is still calculated. This seems to work fine with 1 proc (sequential), but there is some issue in DISTRIBUTE.f where the number of elements on the cap face is set to 0 after partitioning the face among processors (although this doesn't always happen)
Computes weighted surface normal for surface element in deformed configuration. Same as GNNBSURF(), but uses current configuration nodal positions, analogous to GNNBT().
Adding code written to attempt to implement support for a virtual face. This code does not work. I ran into issues with the virtual face being partitioned incorrectly and information being inaccessible when running in parallel. After talking with Vijay, I will instead try to import the face as a member of the boundary condition structure.
Adding code to compute integrals over virtual face. This requires communication between processors when integrals is performed. I have set it up so that master receives all the nodal info it needs to perform the integral, and performs the integral by itself. This could be change so that processors share the work of integration.
In SETBC.f, added a line so that if a face is virtual and there is a Neumann or traction BC defined on it, it will cycle. Otherwise, this leads to a crash.

Also, when testing an idealized LV coupled to a constant flow rate LPN using genBC, I noticed the convergence was suboptimal (less than linear, when it should be quadratic). I think this is because the tangent matrix is inconsistent with the residual vector when a virtual face is used to compute flow rate. In this case, the integrals in the coupled BC tangent matrix contribution (see Moghadam 2013 eq. 27) are not the same. The right one should be taken over the capped surface (endo + virtual cap), while the other should be taken over the uncapped surface). I think this can implemented fairly easily in ADDBCMUL() by cycling in the second loop when the face is virtual. I just need to add and share the virtual flag to lhs%face. Currently, this doesn't compile, but this is what I will fix next.
I believe when a virtual face is used to compute flow rate, the virtual face should have some contribution to the tangent matrix, even if not pressure is applied to it. Because no pressure is applied to it, the tangent matrix contribution should involve an integral over the capped surface and an integral over the uncapped surface. In ADDBCMUL(), I believe this amounts to ignoring the second do loop when the face is virtual. I attempted to do this, but I think the integral (in ValM), as well as the global indices in glob() may be incorrect when a virtual face is used when running in parallel, just like the previous issues I've had to sort out.

Before, the virtual face made no contribution to the tangent matrix because in genBC, I forgot to provide an offset value for the second coupled face (the virtual cap). This meant that the pressure on this face was always zero, so the resistance value for this face was always zero. In FSILS_SOLVE(), if res = 0 for a face, then that face is not considered coupled, so it doesn't contribute anything to the tangent matrix.
Current code produces segfault on Sherlock. Not completely fixed, but see changes.
Reorganizing parallelization code for virtual face handling into functions.

Specifically, in ALLFUN.f, moved communication code for virtual face integration into GatherMasterS and GatherMasterV, which are functions that collect nodal values (like nodal position) onto master, so that master can compute things for a virtual face element, such as normal vector or integration.

Also, in Distribute.f, moved handling for virtual face in PARTFACE() into a function PARTFACEV(), which sets up the parallel data structures for a virtual face.
@LeiShi3374
Copy link
Author

@vvedula22 I pushed a new version. Please check the contact part. Let me know for any questions.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants