Evolving Dirichlet boundary temperature

ASPECT
1

Hello everyone,

Hope you’re doing good. I’m having trouble with boundary temperature model. What I’m doing is mantle convection coupled with core cooling. The core is assumed to be an isothermal heating bath, i.e., we only need to calculate the core temperature (single value, scalar) based on energy balance at core-mantle boundary:
\rho_c C_{p,c} V_c \Delta T = q S \delta t
where \rho_c, C_{p,c}, V_c are core density, specific heat and volume. \Delta T is the temperature change during time interval \delta t. q is heat flux on CMB and S is surface of CMB. All parameters are given or known except heat flux q and \Delta T. Once q is known, temperature change \Delta T is known, and the temperature boundary condition for next step is known, as T_{new} = T_{old} + \Delta T. So basically this is Dirichlet boundary condition but it’s evolving.

How to get q? There are a couple of ways. (1) Adapting dynamic_core.cc where simple conduction heat flux is given in update(), but I realized later that advection also matters at CMB, so dymaic_core may underestimate the heat flux; (2) Making use of Postprocess which prints heat flux at the end of each step; (3) Duplicating functions in Postprocess and adapt them accordingly.

For the 1st way, I have confirmed that the heat flux from my plugin (with adapted dynamic_core inside) is different from that from Postprocess, as expected.

For the 2nd way, by including headers of Postprocess and calling the function aspect::Postprocess::internal::compute_heat_flux_through_boundary_faces(*this), the heat flux can be obtained but it’ still different from Postprocess.

For the last method, funtions from Postprocess are copied to my plugin, including compute_heat_flux_through_boundary_faces() and compute_dirichlet_boundary_heat_flux_solution_vector(). Unfortunately this still gave different heat flux from Postprocess.

There are two more caveats: the heat flux is calculated based on solving (full size?) heat equation on boundary, so this heat flux is accurate, including advection, conduction and more. So I need to get this heat flux. In addition, the Postprocess will be called at the end of each time step, so i have to escape the 1st step (during which initial boundary temperature is known and applied). Starting from 2nd step, I’ll call or use functions from Postprocess in my plugin (my plugin has update() and lives in BoundaryTemperature). I was assuming that the state of the end of n-th step should be the same as that of the beginning of n+1 step. But it turns out that they are different, e.g.,

*** Timestep 561:  t=1.3578e+08 years, dt=36490.2 years
Heat flow at CMB = -29128.7, dT = 0.0118575
   Solving temperature system... 7 iterations.
...
Postprocessing:
     Writing graphical output:           laneuville2013/solution/solution-00458
     RMS, max velocity:                  0.00857 m/year, 0.033 m/year
     Temperature min/avg/max:            250 K, 1607 K, 1969 K
     Heat fluxes through boundary parts: -2.649e+04 W, 4.442e+05 W
     Writing depth average:              laneuville2013/depth_average
     Compositions min/max/mass:          0.07706/1/8.741e+12 // -9.466e-08/3.513e-06/3.943e+05 // -2.527e-08/9.384e-07/1.053e+05 // -2.226e-07/8.259e-06/9.269e+05 // -1.621e-07/6.017e-06/6.754e+05

*** Timestep 562:  t=1.35817e+08 years, dt=36564 years
Heat flow at CMB = -29071, dT = 0.0118579
   Solving temperature system... 7 iterations.
   Solving peridotite system ... 7 iterations.

Note that temperature boundary will be executed before the solve of the system. I have no idea how to handle this difference, any ideas? Deeply appreciate your time and help!

Thanks,
Mingming

2

Hello,
I spent more than a week on this and finally I think I made it. Just declare a static variable to be returned to update() in BoundaryTemperature as heat flow, and a static function to return that variable. Then use heat flow in update().

Thanks,
MM

3

Hi Mingming,

As you have already figured out there are fundamentally different ways to do this, namely:

  1. Use the existing algorithm in postprocess/heat_flux_map.cc
  2. Write a new algorithm in your plugin.

Conceptually reusing code is always preferable, unless there are important reasons not to. The algorithms in postprocess/heat_flux_map have been extensively tested and cover a lot of corner cases, so I would recommend you use them unless they dont work for you.

Concerning the difference in heat flux between postprocessing and the time of update: In the meantime the timestep is advanced and ASPECTs stored solution vectors are shifted by one timestep (and the new solution vector is unknown and left at the same value as the last timestep). Since the CBF method relies on current as well as past timesteps this changes the result and likely leads to a less accurate result at the beginning of a timestep (though how much less accurate I cant estimate, could be small, could be large). I am not sure how this could be prevented. One way would be to store the integrated heat flux during postprocessing (in the postprocessor class). Then access it during `update` from your plugin.

If your implementation with a static variable works, you can certainly use it. For ASPECTs main code we try to avoid static variables, but that is mostly for maintainability. I am sure it can be implemented without static variables, but if it works for you, it works.

Best,

Rene

4

@optimux This is not really about the question you ask, but I’d like to understand this part better as well. Conceptually, the heat flux from the core to the mantle is through an infinitely thin boundary. In truth, there are boundary layers of finite thickness on both sides, but one should understand them to be very thin, and the heat flux is then proportional to the temperature difference \\Delta T between core and mantle outside these boundary layers. But I suspect that these boundary layers are small enough that there is no convection/advection within the boundary layers, just thermal conduction, and that the assumptions of the core_dynamics module are not all that wrong. Would you be willing to elaborate on what you think is wrong from a physical perspective?

Best

Wolfgang

5

Hello Rene,
Thanks for your reply. I gave up using dynamic_core and calling functions in heat_flux_map.cc is that the heat flux difference between Postprocess outputs and the above mentioned ways is getting bigger and bigger, some 20%, you can see the difference above. Although initially the difference is small.

Yes, to be consistent, I’d like to reuse the results from Postprocess and so basically I’m following what you suggested, store the value and then access it in my plugin.

Regards,
MM

6

Hello Wolfgang,

First of all, I don’t think the dynamic_core is wrong, it’s just less accurate. The reason why I did all these is to be consistent with Postprocess.

To answer your question, I do agree that there are boundary layers on both sides of CMB. But from my perspective there are at least 2 occasions where advection might kick in boundary layer. One is the root of plumes, the plumes go up from CMB and thus radial advection should be present. The other could be the case where velocity boundary layer is thinner than thermal boundary layer.
The consistent boundary flux (CBF) used in heat_flux_map.cc ensures that the energy is conserved/balanced rigorously everywhere including the boundaries nearby. Once the temperature field is known, the time derivative term, advection term, source term and diffusion term are known, and so heat flux can be derived following these terms.

That’s what I know. Please be advised that it could be misleading!

Best,
MM

7

I see. Have you thought about comparing the widths of the respective boundary layers? It should be possible to compute their widths based on the material properties.

I cannot say that I have actual knowledge of this, but my intuition tells me that the temperature boundary layer is substantially smaller than the size/diameter of plumes and that heat conduction dominates heat advection. I may of course be totally wrong about that, but I think it would be useful to actually make sure that it’s necessary to fix the boundary condition.

Best
Wolfgang

8

I guess I cannot do this because the core is assumed an isothermal heat bath which has no thermal boundary layer near CMB.

Yes, heat conduction dominates advection in boundary layer. The difference between pure conduction heat flux and CBF one is ~10% (not 20% as I noted before). The size/diameter of plumes could be larger than thermal boundary layer, but the root of the plume (not the head or the neck) might be inside of the thermal boundary layer, though the magnitude of advection would be pretty small, that’s why the difference is not not large.

I’m not an expert in fluid dynamics, so I cannot show you clear examples or analysis :thinking:

Best,
MM

9

Hi everyone,

In previous post, I “shared” a variable among different plugins using static. That was carried out in ASPECT docker image which I believe is built based on Linux platform. When the code is moved to my local installation (macOS), the static way does not apply. Since the plugins are meant to be independent (ASPECT design philosophy), it looks like the plugins cannot talk to each other directly. So I have to created a shared variable (declared as extern) to share the numbers among plugins.

There are 2 more differences between docker image and local installation. One is the way to load symbols during runtime. Docker image(Linux) is more permissible, during runtime, it will load the plugin even with unresolved symbols, and then ld.so will search current plugin, shared libraries, main executable to find the symbol required by this plugin, then link the symbol and execute. However, the linker on macOS, dyld, will not search in main executable. So you will see missing symbol errors during runtime, because there is no hook marked in main executable to link against. That means an independent copy of functions used in Postprocess for heat flux calculation has to be present in my customized plugin.

The other difference is the sensitivity of the code to the mesh. In docker image, I tried many times, for the same parameters, Initial global refinement and Initial adaptive refinement have to be 6 and 3 respectively, while on local installation, it works for 4 and 3 respectively. I’m pretty sure I used same parameters, but why does this happen?

Best,
Mingming

10

@optimux I don’t know whether that will solve your problem, but:

I don’t know what causes your problem, but it might be useful to first try and build through infrastructure that is already there rather than inventing other approaches.

Best
W.

11

Hi @bangerth,

Thanks for the clue, that’s what I’m looking for. I don’t know there is function required_other_postprocessors() I can use directly. I’d like to play with it a little bit.

Best,
MM

12

Hi @bangerth,
Looks like I did not notice at the beginning that your suggestion could work within Postprocess, but I need to ‘share’ or ‘pass’ some numbers from Postprocess to BoundaryTemperature, different namespaces. Is it possible to do so among different namespaces?

Best,
MM

13

Yes. If your boundary temperature object is derived from SimulatorAccess, you can call this->get_postprocess_manager()->template get_matching_postprocessor<...>(). I can’t seem to find a good example, but take a look at the file tests/find_postprocess.cc.

Best
W.

14

That can save a lot, but there’s nothing there. I’ll follow my way first and then switch to this when free.

Thanks,
MM