Soret isothermal wall fix #391

ThomasHowarth · 2024-06-21T16:39:35Z

This PR addresses #389 . Happy to discuss any of the implementation, and let me know if you can see spots for improvement/clarity. Currently, I've not implemented the necessary changes required for EB.

drummerdoc · 2024-06-21T20:30:26Z

Is it possible to check in the test case?

ThomasHowarth · 2024-06-24T08:34:26Z

Source/PeleLMeX_Diffusion.cpp

@@ -917,7 +993,7 @@ PeleLM::differentialDiffusionUpdate(
      GetVecOfConstPtrs(
        getDensityVect(AmrNewTime)), // this triggers proper scaling by density
      GetVecOfConstPtrs(getDiffusivityVect(AmrNewTime)), 0, bcRecSpec,
-      NUM_SPECIES - NUM_IONS, 0, m_dt);
+      NUM_SPECIES - NUM_IONS, 0, m_dt, GetVecOfConstPtrs(spec_boundary));


This probably isn't correct for the boundary MF to be compatible with EFIELD, any suggestions?

Having looked at the mechanisms with ions, the E species has a non-zero Soret coefficient (probably unsurprising), so this probably needs a little bit of thought

I'm open to the idea of just putting an abort for the combination of Efield+Soret+Isothermal and punting on this until someone needs that particular corner case.

Agreed on the punt. :)

baperry2

@ThomasHowarth thanks for figuring this out and adding the test case. After this PR gets merged, I can submit a separate PR to have the case run in testing suite including a check on species balance.

I've been trying out the test case. I got rid of dt_shrink to take larger timesteps and make errors more visible. This fix does drastically reduce, but not quite eliminate, species imbalance errors tracked through temporals/tempSpecies. The remaining imbalance is small enough that I'm not very concerned, but I also don't get why conservation isn't quite exact (it doesn't seem to improve by tightening the diffusion solve tolerances).

A couple comments:

Can the same fix be applied to the EB version of diffuse_scalar? That wouldn't eliminate the issue of isothermal EBs, but would at least allow EB simulations with isothermal domain boundaries and Soret.
It appears that the updated BCs are not applied when diffusion fluxes are computed to compute divu. This should not affect species conservation and with the chi correction on divu, it should not affect the results in the limit of converged SDC iterations.

baperry2 · 2024-06-24T23:05:04Z

Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.H

+  auto eos = pele::physics::PhysicsType::eos();
+  amrex::Real Xt[NUM_SPECIES] = {0.0};
+  amrex::Real a = 0.0;
+  amrex::Print() << prob_parm.fuelID << std::endl;


Remove print statements

baperry2 · 2024-06-24T23:05:37Z

Exec/RegTests/IsothermalSoretChannel/pelelmex_prob.cpp

+  amrex::ParmParse pp_pele("peleLM");
+  std::string fuelName = "";
+  pp_pele.get("fuel_name", fuelName);
+  amrex::Print() << "H2_ID 1 " << H2_ID << std::endl;


remove print statement

ThomasHowarth · 2024-08-07T16:45:07Z

Sorry for the delay, busy month or so...

I've added in the changes so it is computed when computing divu, I've added in the same boundary input to the diffuse_scalar call when running with EB. Code also aborts on Isothermal+Soret+Efield. I am also still seeing some species flux in the species balance files, on the order of 3e-7 or so - any thoughts on this?

baperry2 · 2024-09-04T16:43:56Z

Sorry @ThomasHowarth, I know you've got a few big PRs in the queue waiting for review. I'm going to try to get to them this week.

I've dug into this one a bit and think I figured out the where the remaining conservation error is coming from. With the lagged correction approach, I was only thinking of the implicitly solved $\tilde{D}^{n+1,k+1}$ term in the update equation below, but, as you've now implemented the corrections must also be applied when computing $D^n$ and ${D}^{n+1,k}$, but those are explicitly computed and can't be lagged. Basically, the Soret and Wbar fluxes should be computed first in computeDifferentialDiffusionFluxes and those values should be used to set the inhomogeneous neumann conditions for the diffusive fluxes. For the Soret fluxes that should be easy with some code rearrangement. For the Wbar fluxes it's a little tricky because the Wbar fluxes at the isothermal boundary depend on the species gradients at the boundary, but also affect what the boundary condition needs to be for the species gradients, leading to an implicit definition of the boundary condition. I think it should be okay to compute the Soret terms first, use the boundary species gradients from the Soret terms in computing Wbar fluxes, then use the combined boundary gradients from both Soret and Wbar for the diffusive fluxes. Iterating around the Wbar flux computation would be the fully correct thing though.

Let me know if this makes sense. This is definitely a tricky problem and I'm not sure if I fully understand it. But I did hack a test for the Soret part of the fix (while turning Wbar fluxes off) and that did lead to correct mass conservation.

BTW, I also submitted #410 to make using temporals to track conservation a bit more user friendly.

ThomasHowarth · 2024-09-18T12:50:54Z

Source/PeleLMeX_Diffusion.cpp

@@ -675,36 +784,14 @@ PeleLM::addSoretTerm(
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
    {
-      FArrayBox rhoY_ed;


There was an unnecessary rhoY_ed calculation in the Soret fluxes, which I've now removed

ThomasHowarth · 2024-09-18T13:02:24Z

Source/PeleLMeX_Diffusion.cpp

+        wbarfluxes[lev][idim].setVal(0.0);
+      }
+    }
+    int iter_max = 10;


Currently 10 iterations hardcoded - we could do some sort of residual check instead

10 seems a bit excessive, no?

Probably - based on the test case it converged in 2 or 3 iterations. EDIT: things are blowing up, will keep testing

ThomasHowarth · 2024-09-18T15:55:25Z

Without the Wbar correction, I have also managed to eliminate the conservation error (Up to 10^-23). However, iterating around the wbar fluxes seems to blow up the gradients, and I'm not entirely sure why. I'll keep working on it, but I've also pushed the latest version so you can take a look if you like.

Source/PeleLMeX_Diffusion.cpp

baperry2 · 2024-09-18T17:16:51Z

Nice progress a @ThomasHowarth. This has turned out to be much more involved than expected!

If a solution doesn't jump right out at you for the exploding wbar iterations, do you think we can merge this without the wbar iterations (or maybe hardcoding to 1 "iteration") and then have a separate PR for the final correction? There's a lot here that I don't want to risk falling out of date and even without the last bit it's more correct than what's in there now.

ThomasHowarth · 2024-09-19T14:53:15Z

Source/PeleLMeX_K.H

+  eos.inv_molecular_weight(imw);
+  gradMwmix(i, j, k) = 0.0;
+  for (int n = 0; n < NUM_SPECIES; n++) {
+    imw[n] *= 1000.0; // CGS -> MKS


I thought this was correct, but if I comment out, everything converges nice and quickly - thoughts?

This reflects me being away from things for a while, but the MKS version of mw could be kg/kmole (depending on how it is used elsewhere) so that there is no conversion factor necessary. :)

That is weird because I do think multiplying by 1000 is correct here. When computing MWmix we divide by 1000. Perhaps there's some error cancellation going on?

Thomas Howarth added 11 commits June 21, 2024 15:08

Soret Isothermal wall fix

3c847ae

formatting

145f8e0

formatting

ac393fd

spelling

de385ce

add flag

174348c

fix for efield

f643ab2

formatting

0cd9fab

efield fix

c324e72

ignore eb for now

22d13ff

formatting

0f321fa

bool to int

163dc64

Thomas Howarth and others added 3 commits June 24, 2024 10:27

test case

22a8792

remove autosave file

0b6613f

formatting

7e47b23

ThomasHowarth commented Jun 24, 2024

View reviewed changes

make test case more robust

322a452

baperry2 self-requested a review June 24, 2024 22:26

baperry2 reviewed Jun 25, 2024

View reviewed changes

baperry2 linked an issue Jun 25, 2024 that may be closed by this pull request

Wall boundary conditions (e.g. NoSlipWallIsothermal) do lead to wrong behaviour using the soret effect. #389

Open

Thomas Howarth and others added 10 commits July 2, 2024 18:30

remove prints

f6159d5

remove interaction with efield

d1869c7

fix for divu computation

7a92b9f

formatting

e415fad

remove prints

44f85e3

EB compatibility and Efield fix

db1814b

formatting

236daea

header fix

478467e

missing argument

d439b93

Merge branch 'development' into development

d003712

ThomasHowarth added 2 commits August 12, 2024 20:04

Merge branch 'development' into development

8527f8f

Merge branch 'development' into development

92ba8b0

baperry2 and others added 4 commits September 4, 2024 10:47

Merge branch 'development' into development

fbf6076

Merge branch 'development' into development

70729db

Merge branch 'development' into development

13a8c69

explicit flux calculation, plus wbar fix

3c7d3ed

ThomasHowarth commented Sep 18, 2024

View reviewed changes

removed unnecessary species call in soret term

cea2e63

ThomasHowarth commented Sep 18, 2024

View reviewed changes

Thomas Howarth added 5 commits September 18, 2024 15:09

remove unnecessary captures

fd18c32

unnecessary variable

3178d30

remove commented print

30c0f08

go back to wbar for wbarflux

5011524

formatting

dd5699d

github-advanced-security bot found potential problems Sep 18, 2024

View reviewed changes

Source/PeleLMeX_Diffusion.cpp Fixed Show fixed Hide fixed

ThomasHowarth commented Sep 19, 2024

View reviewed changes

CGS tweak and residual calculation

da1d363

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Soret isothermal wall fix #391

Soret isothermal wall fix #391

ThomasHowarth commented Jun 21, 2024

drummerdoc commented Jun 21, 2024

ThomasHowarth Jun 24, 2024 •

edited

Loading

ThomasHowarth Jun 24, 2024

baperry2 Jun 25, 2024

drummerdoc Jun 25, 2024

baperry2 left a comment

baperry2 Jun 24, 2024

baperry2 Jun 24, 2024

ThomasHowarth commented Aug 7, 2024

baperry2 commented Sep 4, 2024

ThomasHowarth Sep 18, 2024

ThomasHowarth Sep 18, 2024

drummerdoc Sep 18, 2024

ThomasHowarth Sep 18, 2024 •

edited

Loading

ThomasHowarth commented Sep 18, 2024

baperry2 commented Sep 18, 2024

ThomasHowarth Sep 19, 2024

drummerdoc Sep 19, 2024

baperry2 Sep 19, 2024

Soret isothermal wall fix #391

Are you sure you want to change the base?

Soret isothermal wall fix #391

Conversation

ThomasHowarth commented Jun 21, 2024

drummerdoc commented Jun 21, 2024

ThomasHowarth Jun 24, 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

baperry2 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ThomasHowarth commented Aug 7, 2024

baperry2 commented Sep 4, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ThomasHowarth Sep 18, 2024 • edited Loading

Choose a reason for hiding this comment

ThomasHowarth commented Sep 18, 2024

baperry2 commented Sep 18, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ThomasHowarth Jun 24, 2024 •

edited

Loading

ThomasHowarth Sep 18, 2024 •

edited

Loading