change initial condition of RadForce

[chongchonghe]

Description

The RadForce test can be improved to be more realistic and less ad hoc. Instead of initializing with the exact solution and letting the simulation maintain a static state, we start the simulation with a uniform density of rho0 and zero velocity. The left boundary condition is kept the same as before, but the right boundary is set to be foextrap. The end result is exactly the same as before: the radiation blows away the gas and the system comes to the static state given by the differential equation.

Checklist

Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an x inside the square brackets [ ] in the Markdown source below:

• I have added a description (see above).
• I have added a link to any related issues see (see above).
• I have read the Contributing Guide.
• I have added tests for any new physics that this PR adds to the code.
• I have tested this PR on my local computer and all tests pass.
• I have manually triggered the GPU tests with the magic comment /azp run.
• I have requested a reviewer for this PR.

[chongchonghe]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 5 pipeline(s).

[github-actions]

clang-tidy made some suggestions

[chongchonghe]

I made this change in order to do a scaling test for the RHD solver, but the result is surprising: the R1 norm of the error does not change with decreasing dt (see table below). I don't know why this happens.

dt	nx	L1 norm	elasped time
2e8	128	0.0004315252173	4.344377
1e8	128	0.0004296931455	6.812964
5e7	128	0.0004292276377	9.439731
2e8	64	0.0008842796415	2.600089
1e8	64	0.0008833405693	3.727579
2.5e7	64	0.0008808892701	12.810238

[chongchonghe]

I just realized that this problem is in the diffusion limit, so it's not a good test for scaling of the time integration scheme. But in any case, this PR is worth merging.

[BenWibking]

I made this change in order to do a scaling test for the RHD solver, but the result is surprising: the R1 norm of the error does not change with decreasing dt (see table below). I don't know why this happens.

dt nx L1 norm elasped time
2e8 128 0.0004315252173 4.344377
1e8 128 0.0004296931455 6.812964
5e7 128 0.0004292276377 9.439731
2e8 64 0.0008842796415 2.600089
1e8 64 0.0008833405693 3.727579
2.5e7 64 0.0008808892701 12.810238

This will happen if the error is dominated by the spatial discretization rather than the time discretization. In that case, you only get significantly smaller errors if you increase the spatial resolution.

[BenWibking]

For a test problem, I don't think there's anything wrong with initializing to the exact solution. Is there any other benefit to changing this?

[chongchonghe]

For a test problem, I don't think there's anything wrong with initializing to the exact solution. Is there any other benefit to changing this?

I just think a test that is initialized with the exact solution is not super reliable. Support you accidentally changed some of the problem parameters and the simulation finishes at t = 0, or the system is frozen for some reason and never evolves, the test will still pass. This PR shows you can start from a realistic initial condition and end with the exact analytic solution, so why not apply it?

[chongchonghe]

@BenWibking Can you take a look at the setOutflowBoundaryLowOrder boundary condition in test_radiation_force.cpp? Why does it not work? The error message and where it comes from are copied below.

Assertion !S_filled.contains_nan() failed,

	AMREX_ASSERT(!S_filled.contains_nan()); // check ghost zones (usually this is caused by
						// forgetting to fill some components when
						// using custom Dirichlet BCs, e.g., radiation
						// variables in a hydro-only problem)

[BenWibking]

@BenWibking Can you take a look at the setOutflowBoundaryLowOrder boundary condition in test_radiation_force.cpp? Why does it not work? The error message and where it comes from are copied below.

Assertion !S_filled.contains_nan() failed,

	AMREX_ASSERT(!S_filled.contains_nan()); // check ghost zones (usually this is caused by
						// forgetting to fill some components when
						// using custom Dirichlet BCs, e.g., radiation
						// variables in a hydro-only problem)

Ah, I am almost certain the cause is that I wrote this function to only fill the hydro variables. You will have to fill the ghost cells for the radiation variables yourself.

However, I forgot that this was all designed for the full Euler equations, i.e. where there is also an energy equation, which is missing for an isothermal EOS. So your original approach is probably fine. The serious issues only arise when the pressure boundary condition is not specified at an outflow boundary, but here it is not an issue since the pressure is not an independent variable. Somehow I didn't connect the dots here, sorry.

[chongchonghe]

It shouldn't matter anyway, because the gas is definitely outflowing at the upper boundary. In my implementation of setOutflowBoundaryLowOrder, I used something similar to foextrap to pressure (= rho_R * cs^2).

So, the foextrap B.C. is okay? I'll change it back.

[chongchonghe]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 5 pipeline(s).

[BenWibking]

@chongchonghe Can you take a look at the GPU test that failed?

[chongchonghe]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 5 pipeline(s).

[github-actions]

clang-tidy made some suggestions

[chongchonghe]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 5 pipeline(s).

[azure-pipelines]

Azure Pipelines successfully started running 5 pipeline(s).

[github-actions]

clang-tidy made some suggestions

[github-actions]

clang-tidy made some suggestions

[chongchonghe]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 5 pipeline(s).