Boundary_conditions

[kufapetr]

Hello everyone, I implemented my initial and boundary conditions, etc. However, I need to include a change of some low boundary values in density during the calculation. How can I make it? (if there is a possibility to do so). Thanks, Petr

[zingale]

Hi @kufapetr, there's an example of a space- and time-dependent BC in the double_mach_reflection problem (see Exec/hydro_tests/double_mach_reflection.

That problem implements a custom problem_bc_fill.H that fills the ghost cells directly.

[kufapetr]

Hi Michael, thank you for the hint, I'll use it.
Petr

[kufapetr]

Thank you very much (as a late response); this works perfectly. However, I have another similar question: If I need to permanently fill a space- and time-dependent constraint not only to BCs but to either the whole or at least some spatially defined region of the computational domain? What is the best way to do it? Thanks in advance, Petr

On 2023-04-18 15:23, Michael Zingale wrote: Hi @kufapetr [1], there's an example of a space- and time-dependent BC in the double_mach_reflection problem (see Exec/hydro_tests/double_mach_reflection. That problem implements a custom problem_bc_fill.H that fills the ghost cells directly. -- Reply to this email directly, view it on GitHub [2], or unsubscribe [3]. You are receiving this because you were mentioned.Message ID: ***@***.***>

Links: ------ [1] https://github.com/kufapetr [2] #2423 (comment) [3] https://github.com/notifications/unsubscribe-auth/AUJSAEO7SJVSZEXLZGSYY2TXB2IULANCNFSM6AAAAAAXCKZW5U

[maximumcats]

Could you give an example of the kind of constraint you want to apply? Depending on what it looks like, we may need to add some code to support this.

[kufapetr]

Yes, thank you. An example: I need to keep the central region of the computational domain (spherical star) at a constant density, temperature and pressure throughout the entire simulation, regardless of how the surroundings evolve hydrodynamically. Thanks, Petr

On 2023-06-22 13:55, Max Katz wrote: Could you give an example of the kind of constraint you want to apply? Depending on what it looks like, we may need to add some code to support this. -- Reply to this email directly, view it on GitHub [1], or unsubscribe [2]. You are receiving this because you were mentioned.Message ID: ***@***.***>

Links: ------ [1] #2423 (comment) [2] https://github.com/notifications/unsubscribe-auth/AUJSAEKR4Y5F5E3TO5XFZATXMQXDZANCNFSM6AAAAAAXCKZW5U

[maximumcats]

I would recommend achieving this using Prob.cpp and Problem.H in your problem directory to implement the problem_post_timestep() function. This is run after every timestep ends and can be taken advantage of to reset the state to how it was at the beginning of the timestep. Let us know if you need something more strict than that (that is, where the central region would never update at all even in the middle of a timestep).

There are a couple examples of this in Castro/Exec, the one in Exec/science/wdmerger is the most fleshed out. I didn't test this but here is most of a possible implementation I just typed out:

$ cat Problem.H

#ifndef DO_PROBLEM_POST_TIMESTEP
#define DO_PROBLEM_POST_TIMESTEP
#endif

void problem_post_timestep();

$ cat Prob.cpp
#include <Castro.H>
#include <Castro_util.H>

using namespace amrex;

void
Castro::problem_post_timestep()
{
    BL_PROFILE("Castro::problem_post_timestep()");

    MultiFab& S_old = get_old_data(State_Type);
    MultiFab& S_new = get_new_data(State_Type);

    GeometryData geomdata = geom.data();

#ifdef AMREX_USE_OMP
#pragma omp parallel
#endif
    for (MFIter mfi(S_new, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        const Box& bx = mfi.tilebox();

        auto U_old = S_old[mfi].array();
        auto U_new = S_new[mfi].array();

        amrex::ParallelFor(bx,
        [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            Real stellar_radius = 1.0e9_rt;

            GpuArray<Real, 3> r;
            position(i, j, k, geomdata, r);

            Real dist_from_center = 0.0;

            for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
                dist_from_center += (r[dir] - problem::center[dir]);
            }

            dist_from_center = std::sqrt(dist_from_center);

            if (dist_from_center < stellar_radius) {
                for (int n = 0; n < NUM_STATE; ++n) {
                    U_new(i,j,k,n) = U_old(i,j,k,n);
                }
            }
        });
    }
}

[kufapetr]

Hello Max,
thanks for that, it works fine. However, I have one more idea: if I want to reset the profile not as a constant value but as some continuous function (like a sponge function or something) can I use the same principle and then write this function of the old value as the final U_old value?

Thank you very much
Peter

[maximumcats]

In terms of what the software allows, you can make any modification to U_old and U_new that you want in this function, U_old is not read-only. But from a practical point of view I would recommend treating U_old as read-only during a timestep, If you are consistent in how you update U_new, then U_new will become the U_old for the next timestep and you won't need to modify it later.

It is possible I misunderstood your question though, please try asking it again in a different way if so.

[kufapetr]

That is all OK, however, I'm interested in one more thing in this regard: How to get the current TIME here. For example in "problem_bc_fill.H" the current time is loaded, I tried to get it here but so far I can't.

[kufapetr]

Dear all, I apologize I ask once again for the same thing, but I accidentally deleted all the previous incoming messages which I did not backup; could you please send me again the supporting code for this (or an example reference)... Thank you very very much, Petr

On 2023-06-23 04:52, Petr_Kurfurst wrote: Yes, thank you. An example: I need to keep the central region of the computational domain (spherical star) at a constant density, temperature and pressure throughout the entire simulation, regardless of how the surroundings evolve hydrodynamically. Thanks, Petr On 2023-06-22 13:55, Max Katz wrote: > Could you give an example of the kind of constraint you want to apply? > Depending on what it looks like, we may need to add some code to > support this. > > -- > Reply to this email directly, view it on GitHub [1], or unsubscribe > [2]. > You are receiving this because you were mentioned.Message ID: > ***@***.***>

Links: ------ [1] #2423 (comment) [2] https://github.com/notifications/unsubscribe-auth/AUJSAEKR4Y5F5E3TO5XFZATXMQXDZANCNFSM6AAAAAAXCKZW5U

[zingale]

I'm not sure what examples you mean, but the previous examples are still in this discussion thread on github.

[kufapetr]

Oh, sure, I have it. Sorry for bothering you. Petr

[kufapetr]

Hello everyone, I have another question: In the inputs file, I limit the minimum density for example as "castro.small_dens = 1.e-18" and at the same time in the inputs file, I also specify "castro.retry_small_density_cutoff" which according to the manual should prevent retry in low-density areas. However, after some time the negative density message will still appear, along with "Advance was unsuccessful with reason: invalid density; proceeding to a retry". What's the problem? Should the limiting commands be entered somewhere other than in the "inputs" file? Thank you, Petr.

Translated with www.DeepL.com/Translator (free version)

[kufapetr]

Or, as I see in some examples, there must be a specified value at the "castro.retry_small_density_cutoff", I will run it once again! Petr

[zingale]

Hi Petr, indeed, this is discussed a bit here: https://amrex-astro.github.io/Castro/docs/timestepping.html#index-1

The reason why density could be going negative is that you are not using enough resolution, so some steep gradient is unresolved. If you don't care about this feature, then you could set castro.retry_small_density_cutoff to be a value larger than castro.small_dens that is still small enough that you don't care about it.

Another approach is to try the algorithm that adjusts the fluxes to try to prevent hydro from producing negative densities. See this:

https://amrex-astro.github.io/Castro/docs/Hydrodynamics.html#flux-limiting

[kufapetr]

Hello, everybody,
under certain initial conditions, the calculation breaks down after a certain time for the appearance of NaNs in rho_E. What are the effective ways to fix this problem; note that I have set values for castro.small_dens, castro.small_pres, castro.small_temp, eos.eos_input_is_constant = 0. Is there any possibility, for example, to control the values of rho_E directly?

[kufapetr]

one more thing: I lowered the cfl number to 0.4 and so far it seems to have helped.

[maximumcats]

In principle you can do whatever you want, we have a hook that allows you to modify the state at the end of each timestep. However, we also shouldn't be letting you get into a state where you are generating NaN's. Can you compile in debug mode, and use amrex.fpe_trap_invalid=1 (more info on our debugging page), and let us know what the content of the Backtrace file is when you crash?

[kufapetr]

Dear Max, thank you for the reply, I'll do the debugging and let you know. However, this will take some time because the calculation in the debug mode is significantly slower and the crash may occur in a rather late phase of computation. Best, Petr

On 2023-10-02 14:28, Max Katz wrote: In principle you can do whatever you want, we have a hook that allows you to modify the state at the end of each timestep. However, we also shouldn't be letting you get into a state where you are generating NaN's. Can you compile in debug mode, and use amrex.fpe_trap_invalid=1 (more info on our debugging page [1]), and let us know what the content of the Backtrace file is when you crash? -- Reply to this email directly, view it on GitHub [2], or unsubscribe [3]. You are receiving this because you commented.Message ID: ***@***.***>

Links: ------ [1] https://amrex-astro.github.io/Castro/docs/debugging.html [2] #2423 (reply in thread) [3] https://github.com/notifications/unsubscribe-auth/AUJSAEIRTC4ZGYIYJG6R3ILX5KXOTAVCNFSM6AAAAAAXCKZW5WVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TCNRUGUYTK

[zingale]

run with the current executable but increase the checkpoint frequency to be something like every 10 steps or so. Then once it crashes, rebuild in debug mode and restart from the last checkpoint file, and you should only need to evolve for (at most) 10 steps to get the crash.

[kufapetr]

Hello everybody,
I'm planning radiative-hydrodynamic simulations of a supernova explosion, focusing on the supernova's interactions with its surroundings. Which of the Castro templates is best to use as a starting point, which one is closest to this intention? I don't see a supernova explosion directly there. Petr

[zingale]

what type of supernova? Not that neutrino radiation was removed some time ago from the code.

[kufapetr]

Standard type II, let's say RSG 15M_sun, 1000R_sun. I can use the MESA model for the initial state as an input file and, if necessary, I can also use the very early phase of SN blast calculation (until the shock breakout) from 1D SNEC pre-calculation and remap it into Castro. Neutrino radiation is not necessary. Thanks, Petr Initial profile On 2023-10-09 14:43, Michael Zingale wrote:

what type of supernova? Not that neutrino radiation was removed some time ago from the code. -- Reply to this email directly, view it on GitHub [1], or unsubscribe [2]. You are receiving this because you commented.Message ID: ***@***.***>

Links: ------ [1] #2423 (comment) [2] https://github.com/notifications/unsubscribe-auth/AUJSAEMTIW2DSTMOTNBEY6LX6PWQDAVCNFSM6AAAAAAXCKZW5WVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TEMZQGEZDG

[zingale]

I am not sure if we have a setup that is a match for what you want to do. But you can start with Castro/Exec/science/subchandra. That puts a white dwarf on the grid and uses full Poisson gravity. You might want to change to a multipole. It has reactions enabled, but not radiation.

[kufapetr]

Hello all, I understand how to interpolate a 1D model using Model Parser, even on the part of the computational domain; it is OK and works perfectly. However, I cannot figure out how to proceed if I want to remap the initial state of some quantities, like density, temperature, etc, from a multidimensional (like 2D) precomputed model stored in some "model file"? It doesn't, in fact, have to be an interpolation; it can just be a mapping of the values onto the Castro grid within an identical grid geometry.

[zingale]

we don't have an example of that capability right now, so I think that you would have to write your own reader and interpolator/mapper.

[kufapetr]

Of course, I can write my own reader/mapper but where it should be located? For example, the model parser is located in Util/model_parser_cxx/ where in file model_parser.H (if I understand it correctly) is the interpolator. It is, however, linked elsewhere - probably to AMREX core, etc. Where should I place the mapper to be called correctly?

[zingale]

it's probably easiest that you start just by doing it in your problem directory. For each file that needs to be compiled, you'll have to add a line in the Make.package in that directory.

[kufapetr]

Hello, I'm interested in one more thing: How to get the current TIME in Prob.cpp and Problem.H in my problem directory, where I implement the problem_post_timestep() function; I need to achieve this function being on or off in various time intervals. For example in "problem_bc_fill.H" the current time is loaded, I tried to get it in Prob.cpp but so far I can't. Thanks a lot.

[zingale]

In Prob.cpp, you can get the time as:

Real time = state[State_Type].curTime();

see the version in Exec/science/wdmerger for an example.

[kufapetr]

Hello, I have the following question: I calculate a standard hydro model (in 2D) and I need to integrate (not permanently but at least at certain selected times) the total mass inside a region, namely inside a selected radius. So, not within the whole domain but only within a certain part of it. What is the best way to do it please?

[zingale]

hi, can you start a new discussion thread for each different topic -- this one still refers to boundary conditions, so it is a bit confusing to follow.

[zhuohongbin]

Hello, I also want to implemented inflow boundary conditions, etc. For example, constant heat flux in low boundary. How can I make it
or if there is any example to do so. Thanks, ZHUO

[zingale]

you can take a look at Castro/Exec/science/Detonation. There is a problem_bc_fill.H file there that fills the ghost cells with an ambient thermodynamic state. That should be a good starting point for you.