Fixed to the phase-field model #515

shkodm · 2024-05-27T23:29:20Z

Synchronising various fixes I did to the phase-field mode over the past months:

Split initialisation of wall normals in 2 parts, so it is possible to supply user defined wall normals using RunR
Fix various edge cases and bugs with wetting boundary conditions
Fix pressure boundary conditions
Add markers for the saturation

Pressure boundary conditions fixes:

Now use proper ZouHe boundary conditions, with 3 variants. (Some can be removed later)
- fflux enforces constant flux -> Breaks when interface enters the outlet - default (i.e with WPressure or EPRessure name)
- fixed phase-field (e.g. EPressurefpf) makes sure the inlet / outlet phase-field is the same as specified (i.e. sum of distributions forced to be equal to it). Effect: interface propagation stop at the outlet, "impossible-to-displace liquid". But seems to not work with the inlet
- open suffix (i.e. EPressureopen) does not enforce the phase-field. Interface propagates into the outlet, splits in 2,
  there are so strange effects, but otherwise does not look too crazy
Gradient / Laplacian on the velocity / pressure boundaries is always zero (does not touch the other side of the domain)
Wall normal calculations do not take into account the other side of the domain (periodicity is switched off). This works only for the case of X direction and when the flag InvasionDrainage is non-zero. And when the boundary conditions are the first / last layers of the domain. Otherwise, one has to compile the model with autosym option, and mark the inlet / oulets with Symmetry boundaries, so that the periodicity is switched off
InvasionDrainage marks if the case is drainage (non-wetting displacing wetting), or invasion (wetting displacing non-wetting). This sets some densities and phase-field on the inlet / outlet
Apply boundary conditions before VTK export, so the inlet / outlet zones don't look strange due to streaming

Merging with master, because I (and I think other people who are using phase-field model are mostly using master branch and the fixes are quite crucial)

Now doing some testing to make sure I did not mess up when refactoring. So draft PR for now.
cc @TravisMitchell

Compute pf_f in the normal direction as soon as possible, because we probably want to use it later in the computations

This way we can both initialise external and internally based on some approximations

TravisMitchell · 2024-05-27T23:39:58Z

models/multiphase/d3q27_pf_velocity/Boundary.c.Rt

+<?R
+for (ii in 1:cube_faces){
+for (j in 1:length(pressure_boundary_types)) {  ?>
+CudaDeviceFunction void <?%s paste0(C(my_pressure_boundaries[ii]), pressure_boundary_types[j]) ?>(){
    real_t d = getRho();


@shkodm should d be based off a prescribed phase on the boundary rather than from getRho - which I believe calls from PhaseF field

yea, it is then set properly in the if condition below

I did see this, but figured it was maybe not best to define it with a value knowing it will be overwritten. I was also curious how this worked with the open condition - but perhaps I need to just look at the code a bit deeper. I will try to do this before tomorrow when we meet but have a bit of a marking pile with end of semester.

llaniewski · 2024-05-28T01:02:19Z

models/multiphase/d3q27_pf_velocity/Boundary.c.Rt

 }

 CudaDeviceFunction real_t getSpecialBoundaryPoint() {
 	return IsSpecialBoundaryPoint;
 }

+/* Calculate the actual wall norm based on the solid
+ *field */
+CudaDeviceFunction real_t Init_real_wallNorm()


This function should be void (shouldn't be returning anything).

llaniewski · 2024-05-28T01:37:01Z

models/multiphase/d3q27_pf_velocity/Boundary.c.Rt

+<?R
+for (ii in 1:cube_faces){
+for (j in 1:length(pressure_boundary_types)) {  ?>
+CudaDeviceFunction void <?%s paste0(C(my_pressure_boundaries[ii]), pressure_boundary_types[j]) ?>(){


remove the "C(", as it doesn't make sense here.

llaniewski · 2024-05-28T01:46:57Z

models/multiphase/d3q27_pf_velocity/Boundary.c.Rt


-CudaDeviceFunction void <?R C(my_pressure_boundaries[ii])?>(){
+<?R
+for (ii in 1:cube_faces){


Could we change this to a pattern of form table with expand.grid, add relevant stuff as columns and use this table everywhere? It is in other places of code and I think is more clear.

BCs = expand.grid(side=1:4, type=c("Pressure","Velocity"), subtype=c("fix_flux","fix_value")) BCs$side_letter = c("W","E","S","N")[BCs$side] BCs$name = paste0(BCs$side_letter, BCs$type, "_", BCs$subtype)

and then in Dynamics.R one can do:

AddNodeType(BCs$name, group="BOUNDARY")

and in Boundary.c.Rt:

<?R for (bc in rows(BCs)) { ?> CudaDeviceFunction void <?%s bc$name ?>() { <?R ... n=my_normal_directions[bc$side,] ... ?> } <?R } ?>

codecov · 2024-05-28T03:02:54Z

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 36.36%. Comparing base (2362a9d) to head (98accf4).
Report is 24 commits behind head on master.

❗ Current head 98accf4 differs from pull request most recent head 9a0f8d6

Please upload reports for the commit 9a0f8d6 to get more accurate results.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #515      +/-   ##
==========================================
- Coverage   44.24%   36.36%   -7.89%     
==========================================
  Files         166      159       -7     
  Lines        7567     6790     -777     
==========================================
- Hits         3348     2469     -879     
- Misses       4219     4321     +102

Flag	Coverage Δ
d2q9	`?`
d3q27_PSM_NEBB	`36.36% <ø> (+0.01%)`	⬆️
d3q27_pf_velocity	`?`

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

shkodm · 2024-06-04T23:13:13Z

@nparaskevas you can test it now if you want. I provided the description of different boundary conditions variants in the PR description. Just using WPressure on the inlet and EPressurefflux or EPressurefpf or EPressureopen on the outlet would work, but you should play around.

Don't forget to add InvasionDrainage value="2" to the setting in case of drainage, or value = 1 in case of invasion

Need to normalise by the sum of weights obviously (otherwise would leads to crazy low values of the phase field in random points of the domain)

TravisMitchell · 2024-09-05T00:46:22Z

models/multiphase/d3q27_pf_velocity/Boundary.c.Rt

+    if ((NodeType & NODE_BOUNDARY) == <?%s paste0("NODE_EPressure", pressure_boundary_types[j]) ?>){
+       d = InvasionDrainage == 2? Density_h : Density_l;
+       myPhaseField = InvasionDrainage == 2? 1.0 : 0.0;
+    }



@shkodm I think here you should be setting:
PhaseF = PhaseField; d = getRho();
You then shouldn't need myPhaseField as a variable in this function? I'm not sure what your InvasionDrainage setting here is checking as you specify it should be 1 or -1 in the Dynamics.R. This is causing an issue (i think) for pressure-pressure cases where fpf is used as it effectively forces PhaseField=0 on all pressure boundaries.

@TravisMitchell The documentation in Dynamics.R is incorrect (copypasted), InvasionDrainage should have values 1 or 2. d should be set explicitly because I don't want any streamed values obtained from getRho() This was the whole point of those if statements, without it would not work properly.

@shkodm, wouldn't it be much simpler/neater to take PhaseField (zonal setting defined in xml) and then interpolate for d based off Density_h and Density_l. This would mean that the myPhaseField would be simply set from the zonal condition applied in the xml rather than relying on this InvasionDrainage setting.

@TravisMitchell yea, interpolation is a good idea. I think I implemented it with an additional setting originally, because InvasionDrainage (or whatever the name of the flag) is still necessary to know if we want to "reflect/mirror" the boundaries on the inlet (to have proper wall normal not pointing to the other side of domain), or if we want to keep periodic boundaries.

Anyway, I will remove the if statements here later today or tomorrow and double check the boundary conditions.

cc @nparaskevas

Dmytro Sashko and others added 12 commits January 16, 2024 08:46

Add fixes for building on macOS

370a3c7

Adjust boundary computation for angle 90 (or other similar cases)

c3ac8f4

Compute pf_f in the normal direction as soon as possible, because we probably want to use it later in the computations

Split normal initialisation in two part

c3b563f

This way we can both initialise external and internally based on some approximations

Fix geometric model

5bcc8d4

Merge remote-tracking branch 'origin/master' into fractures_master

1d9a6f3

Fix undified behaviours

fd057e6

Add utilities

1b4eca6

Fix weighting

86d7854

WIP Pressure boundary condition fix

731ac87

Add Lukaszes boundary conditions

c45c997

Rewrite fixes generically

e20e665

Add comments

fe4f08b

TravisMitchell reviewed May 27, 2024

View reviewed changes

llaniewski requested changes May 28, 2024

View reviewed changes

Make more explicit

8c48c8d

shkodm added 3 commits May 28, 2024 13:05

Fix

34ea425

Fix

98accf4

Change the default type

660ce39

shkodm added 4 commits June 14, 2024 11:19

Add setting pressure on the boundaries using RunR

b4696c3

Fix typo

9a0f8d6

Fix normalization for the special boundary points

5ca385b

Need to normalise by the sum of weights obviously (otherwise would leads to crazy low values of the phase field in random points of the domain)

Add handling of errors arising in RunPython

33a2b7f

TravisMitchell reviewed Sep 5, 2024

View reviewed changes

shkodm added 3 commits September 5, 2024 15:54

Improve the phase field and density on Boundaries

d69ce1f

Fix typo

b34c5b4

Add more robust version of special boundary points

9d9e583

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fixed to the phase-field model #515

Fixed to the phase-field model #515

shkodm commented May 27, 2024 •

edited

Loading

TravisMitchell May 27, 2024

shkodm May 27, 2024

TravisMitchell May 27, 2024

llaniewski May 28, 2024

llaniewski May 28, 2024

llaniewski May 28, 2024

codecov bot commented May 28, 2024 •

edited

Loading

shkodm commented Jun 4, 2024

TravisMitchell Sep 5, 2024 •

edited

Loading

shkodm Sep 5, 2024

TravisMitchell Sep 5, 2024

shkodm Sep 5, 2024 •

edited

Loading

Fixed to the phase-field model #515

Are you sure you want to change the base?

Fixed to the phase-field model #515

Conversation

shkodm commented May 27, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

codecov bot commented May 28, 2024 • edited Loading

Codecov Report

shkodm commented Jun 4, 2024

TravisMitchell Sep 5, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

shkodm Sep 5, 2024 • edited Loading

Choose a reason for hiding this comment

shkodm commented May 27, 2024 •

edited

Loading

codecov bot commented May 28, 2024 •

edited

Loading

TravisMitchell Sep 5, 2024 •

edited

Loading

shkodm Sep 5, 2024 •

edited

Loading