Velocity boundary condition issues in a 3D model with 'box with lithosphere' model

[maziqi96]

Dear developers,
I hope you’re doing well. I have successfully set up a 2D model with a ‘box with lithosphere’ model and am working on extending the same boundary conditions to a 3D box model shown below.

My velocity boundaries setup are shown below:

subsection Boundary velocity model
set Prescribed velocity boundary indicators = 4 z: function, 6 x: function, 7 x: function
  subsection Function
    set Variable names          = x,y,z
    set Function constants      = vx_top=0.05, vx_bot=0.0, vy_top=0, vy_bot=0; vz_bot=0.0, vz_top=0.0, x0=1320e3,y0=1320e3, r_p=200e3, rho_m=3400, alpha=3e-5, g=9.8, delta_T=250, eta0=1e21, readj=11.4524e-4,gamma=3.1
    set Function expression     = if(z==0e3, vx_bot, vx_top);0;if(z==0e3,if(((x-x0)^2+(y-y0)^2)<r_p^2,0.18,0),vz_top)
  end
end

However, errors occurred when I tried to run the model above:

ERROR: FunctionParser failed to parse
	'Boundary velocity model.Function'
with expression
	'if(z==0e3, vx_bot, vx_top);0;if(z==0e3,if(((x-x0)^2+(y-y0)^2)<r_p^2,0.18,0),vz_top)'More information about the cause of the parse error 
is shown below.
An error occurred in line <128> of file </home/mazq/software/aspect/tmp/unpack/deal.II-v9.4.0/source/base/parsed_function.cc> in function
    void dealii::Functions::ParsedFunction<dim>::parse_parameters(dealii::ParameterHandler&) [with int dim = 3]
The violated condition was: 
    this_c.size() == 2
Additional information: 
    Invalid format

Stacktrace:
-----------
#0  /home/mazq/software/aspect/deal.II-v9.4.0/lib/libdeal_II.g.so.9.4.0: dealii::Functions::ParsedFunction<3>::parse_parameters(dealii::ParameterHandler&)
#1  /home/mazq/software/aspect/aspect/build/aspect: aspect::BoundaryVelocity::Function<3>::parse_parameters(dealii::ParameterHandler&)
#2  /home/mazq/software/aspect/aspect/build/aspect: aspect::BoundaryVelocity::Manager<3>::parse_parameters(dealii::ParameterHandler&)
#3  /home/mazq/software/aspect/aspect/build/aspect: aspect::Simulator<3>::Simulator(ompi_communicator_t*, dealii::ParameterHandler&)
#4  /home/mazq/software/aspect/aspect/build/aspect: void run_simulator<3>(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, bool, bool, bool)
#5  /home/mazq/software/aspect/aspect/build/aspect: main
--------------------------------------------------------

Aborting!
----------------------------------------------------

Could you please help me check where went wrong? I would greatly appreciate it!

Cheers,
Ziqi

[jbnaliboff]

Hi Maziqi,

After reformatting your function velocity statement, I also don’t see exactly what the issue. I think what you have should work if the model is indeed setup as 3D.

if(z==0e3, vx_bot, vx_top); \
0; \
if ( z==0e3, \
     if( ( (x-x0)^2 + (y-y0)^2 ) < r_p^2, 0.18, 0), \ 
     vz_top)

Maybe others can identify what the issue is? To help diagnose, can you also post the mesh and geometry sections of the PRM files, and confirm that the model is set to be 3D and not 2D?

Cheers,
John

[bangerth]

@maziqi96
The error is in your list of constants:

vx_top=0.05, vx_bot=0.0, vy_top=0, vy_bot=0; vz_bot=0.0, vz_top=0.0, x0=1320e3,y0=1320e3, r_p=200e3, rho_m=3400, alpha=3e-5, g=9.8, delta_T=250, eta0=1e2

This should be a list of name=value separated by commas, but you have a semicolon in one place.

Best
W.

[maziqi96]

Hi Wolfgang and John,
Thank you very much for your timely and helpful reply! You’re right! Well-spotted! My model works after correcting that rookie mistake. Thank you!

Cheers,
Ziqi

[maziqi96]

Hi,
Thanks again for helping me spot the typo earlier. However, the model couldn’t converge after the first timestep. Model results show that the pressure and viscosity are not what we want them to be.

Below are the 2D pressure and viscosity snapshots, which we want to extend to our 3D model.

I have attached the prm file and the slurm file below. Could you please help me check what’s wrong with my model setup? I would greatly appreciate it!
original.prm (8.3 KB)
slurm-5855402.txt (8.8 KB)

Cheers,
ziqi

[jbnaliboff]

@maziqi96

Could you please help me check what’s wrong with my model setup? I would greatly appreciate it!

I’m happy to help with debugging the issue, but to start it would be very helpful to know what you expect the output (viscosity, temperature, etc) to look like.

A few other recommendations to help with debugging:

1. Start with a 2D only model, and then expand it to 3D after you are confident your model setup is correct in 2D.
2. Once you move to 3D, begin with a very narrow 3D model so that it effectively reproduces the same (or very similar) results as your 2D models.
3. To help debug issues with the initial composition and temperature setup, try using the no Advection, no Stokes nonlinear solver scheme. This solver scheme will only setup the initial conditions and run the post processors, which allows quickly checking if your setup looks correct to first order. For the visco plastic material model, the Reference strain rate parameter is used for calculating the dislocation creep and plastic viscosities.

Cheers,
John

[maziqi96]

Hi John,
Thank you very much for your reply! We do have a working 2d model. I will try the second and third suggestions. Thank you!

Cheers,
Ziqi

[maziqi96]

Hi John,
I hope this letter finds you well. I want to thank you again for your previous helpful suggestions. I have been testing 3d models according to your suggestions. I did get working models after I narrowed down the width of the 3d model (original.prm (8.2 KB)). However, it only worked for 27 Myrs and then crashed due to the following error:

 *** Timestep 385:  t=2.756e+07 years, dt=71590.9 years
   Solving mesh surface diffusion
   Solving mesh displacement system... 10 iterations.
   Solving temperature system... 15 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 28+0 iterations.

   Postprocessing:
     Temperature min/avg/max: -1.633 K, 1812 K, 2432 K 
     Topography min/max:      -409.9 m, 0.7426 m



----------------------------------------------------
Exception on MPI process <30> while running mesh refinement plugin <N6aspect14MeshRefinement9ViscosityILi3EEE>: 

--------------------------------------------------------
An error occurred in line <970> of file </work/n03/n03/mazq/modules/aspect/source/material_model/utilities.cc> in function 
    double aspect::MaterialModel::MaterialUtilities::average_value(const std::vector<double, std::allocator<double> >&, const std::vector<double, std::allocator<double> >&, const aspect::MaterialModel::MaterialUtilities::CompositionalAveragingOperation&)
The violated condition was: 
    parameter_values[i] > 0
Additional information: 
    All parameter values must be greater than 0 for harmonic averaging!

Stacktrace:
-----------
#0  /work/n03/n03/mazq/modules/aspect/release/aspect: aspect::MaterialModel::MaterialUtilities::average_value(std::vector<double, std::allocator<double> > const&, std::vector<double, std::allocator<double> > const&, aspect::MaterialModel::MaterialUtilities::CompositionalAveragingOperation const&)
#1  /work/n03/n03/mazq/modules/aspect/release/aspect: aspect::MaterialModel::ViscoPlastic<3>::evaluate(aspect::MaterialModel::MaterialModelInputs<3> const&, aspect::MaterialModel::MaterialModelOutputs<3>&) const
#2  /work/n03/n03/mazq/modules/aspect/release/aspect: aspect::MeshRefinement::Viscosity<3>::execute(dealii::Vector<float>&) const 
#3  /work/n03/n03/mazq/modules/aspect/release/aspect: aspect::MeshRefinement::Manager<3>::execute(dealii::Vector<float>&) const
#4  /work/n03/n03/mazq/modules/aspect/release/aspect: aspect::Simulator<3>::refine_mesh(unsigned int)
#5  /work/n03/n03/mazq/modules/aspect/release/aspect: aspect::Simulator<3>::run()
#6  /work/n03/n03/mazq/modules/aspect/release/aspect: void run_simulator<3>(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, bool, bool, bool)
#7  /work/n03/n03/mazq/modules/aspect/release/aspect: main
--------------------------------------------------------

Then I tried using the ‘no Advection, no Stokes’ nonlinear solver scheme, which did help prevent the model from crashing but the velocity field doesn’t work the way we want it to.
What we want the bottom velocity and bottom temperature to be is like:
Bottom_inflow.pdf (75.5 KB)

But in the ‘no Advection, no Stokes’ model, where everything else is the same:

diff test1.prm ../compare_2D_3D/test9.prm 
11c11
< set Output directory                       = output_1
---
> set Output directory                       = output_9
17d16
< set Nonlinear solver scheme 		   = no Advection, no Stokes

The velocity field is 0 everywhere in the ‘no Advection, no Stokes’ model like this:

Could you please help me check where could go wrong to cause this when you’re available? That will be much appreciated. Thank you!

Cheers,
Ziqi

[bangerth]

@maziqi96 What the error message says is that if you are “averaging” parameters in the material model using the “harmonic averaging” method, then the values must all be positive. This makes sense given the formula used for averaging. (If you’re unclear what “averaging” means, take a look at Section 3.3 of https://www.math.colostate.edu/~bangerth/publications/2017-boussinesq.pdf, or the manual). I don’t know what material parameter it is that becomes negative (could be the density, could be the viscosity); your next task is to find out why it is happening, perhaps by running the program in a debugger.

As for why “no Stokes, no advection” does not work: If you select this solver, you are simply not solving for anything at all. The solution will simply never change. This is of course not useful for actual applications, but sometimes useful for debugging other kinds of problems.

Best
W.

[maziqi96]

Hi Wolfgang,
Thanks for your timely and helpful response. I got it! I have read the paper about averaging and ran the model under debug mode to figure out which parameter becomes negative. Thank you very much!

Cheers,
Ziqi

[maziqi96]

Hi Wolfgang,
After I ran it in debug mode, I found that it is the diffusion creep (temperature) value becomes negative.

An error occurred in line <100> of file </work/n03/n03/mazq/modules/aspect/source/material_model/rheology/diffusion_creep.cc> in function
    double aspect::MaterialModel::Rheology::DiffusionCreep<dim>::compute_viscosity(double, double, unsigned int, const std::vector<double, std::allocator<double> >&, const std::vector<unsigned int>&) const [with int dim = 3]
The violated condition was:
    viscosity_diffusion > 0.0
Additional information:
    Negative diffusion viscosity detected. This is unphysical and should
    not happen. Check for negative parameters. Temperature and pressure
    are -1.6328450332290574 K, 522753992.36101103 Pa.

I would appreciate it if you could give me some suggestions on what might cause this. Thank you very much!

Cheers,
Ziqi

[bangerth]

Ziqi,
I’m afraid I will have to leave this to others who may know more about the material models.

That said, one message is already clear: Always first run in debug mode

Best
W.

[jbnaliboff]

Hi Ziqi,

I’m surprised running in release mode did not produce the same error, but that is a separate issue.

Can you check via Paraview/Visit or loading your data into Python to where the temperature is negative?

Likewise, can you confirm that the initial temperature distribution is what you expect it to be? Looking at your initial temperature function from the prm file, it appears there are areas where the temperature might be zero?:

  subsection Function
    set Variable names      = x,y,z
    set Function constants  = DeltaT=250, r_p=200e3,x0=1000e3,y0=0e3
    set Function expression = if(z==0,if(((x-x0)^2+(y-y0)^2)<r_p^2, DeltaT*(1-((x-x0)^2+(y-y0)^2)/r_p^2),0),0)
  end

Is this indeed the case (initial temperature values equal to 0)?

However, in the boundary temperature model you have

subsection Boundary temperature model
  set Fixed temperature boundary indicators = 0,1,2,3,4,5
  set List of model names = initial temperature

  subsection Initial temperature
    set Maximal temperature = 2300
    set Minimal temperature = 273
  end
end

set Minimal temperature = 273 at first glance seems inconsistent with your initial temperature function, and I’m not sure what this combination of parameters will produce.

There are always small over/undershoots that develop when solving the energy equation, and one way to make sure your temperature values are never below 0 is to have boundary temperatures and initial temperatures no lower than 273 (K). It’s possible I’m misreading what the initial temperature function should produce, but I would also be surprised if you have a propagating temperature undershoot error > 273 K.

Long story short, I think you likely need to revise your initial and boundary temperature specifications.

Cheers,
John

[bangerth]

The assertion in debug mode terminates the program before it comes to the
place where it crashes in release mode.

Best
W.