Skip to content

[WIP] new melt formulation #6344

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 7 commits into
base: main
Choose a base branch
from
Open

[WIP] new melt formulation #6344

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
87bdd52
first try new melt formulation
jdannberg Jun 11, 2025
135e76f
make regularization an input
jdannberg Jun 11, 2025
cb22c53
fix preconditioner
jdannberg Jun 13, 2025
6ae67ea
remove is melt cell and add a test
jdannberg Jun 13, 2025
207255d
pressure_normalization_and_do_not_scale_regularization_term
jdannberg Jun 18, 2025
db56e0b
add arbogast benchmark plotting
jdannberg Jun 19, 2025
0d0a113
use inverse compaction viscosity
jdannberg Jun 19, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 2 additions & 22 deletions include/aspect/melt.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -219,28 +219,6 @@ namespace aspect
* permeability divided by fluid viscosity. Units: m^2/Pa/s.
*/
virtual double reference_darcy_coefficient () const = 0;

/**
* Returns the cell-averaged and cut-off value of p_c_scale,
* the factor we use to rescale the compaction pressure and to
* decide if a cell is a melt cell.
* The last input argument @p consider_is_melt_cell determines if
* this computation takes into account if a cell is a "melt cell".
* Melt cells are cells where we solve the melt transport equations,
* as indicated by the entries stored in the is_melt_cell vector of
* the melt handler. In case @p consider_is_melt_cell is set to true,
* this function returns a value of zero if the cell is not a melt cell.
* If @p consider_is_melt_cell is set to false the computation
* disregards the information about which cells are melt cells,
* and computes p_c_scale from the cell-averaged Darcy coefficient
* for all cells. This is needed for example when we want to update
* the is_melt_cell vector and need to find out which cells should be
* marked as melt cells.
*/
double p_c_scale (const MaterialModel::MaterialModelInputs<dim> &inputs,
const MaterialModel::MaterialModelOutputs<dim> &outputs,
const MeltHandler<dim> &melt_handler,
const bool consider_is_melt_cell) const;
};


Expand Down Expand Up @@ -416,6 +394,8 @@ namespace aspect
* be averaged cell-wise.
*/
bool average_melt_velocity;

double regularization;
};
}

Expand Down
23 changes: 8 additions & 15 deletions source/postprocess/visualization/melt.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,9 +59,9 @@ namespace aspect
if (property_name == "fluid density gradient")
for (unsigned int i=0; i<dim; ++i)
solution_names.emplace_back("fluid_density_gradient");
else if (property_name == "compaction pressure")
else if (property_name == "total pressure")
{
solution_names.emplace_back("p_c");
solution_names.emplace_back("p_t");
}
else
{
Expand Down Expand Up @@ -130,11 +130,6 @@ namespace aspect
AssertThrow(melt_outputs != nullptr,
ExcMessage("Need MeltOutputs from the material model for computing the melt properties."));

const double p_c_scale = Plugins::get_plugin_as_type<const MaterialModel::MeltInterface<dim>>(this->get_material_model()).p_c_scale(in,
out,
this->get_melt_handler(),
true);

for (unsigned int q=0; q<n_quadrature_points; ++q)
{
unsigned output_index = 0;
Expand All @@ -156,16 +151,14 @@ namespace aspect
}
--output_index;
}
else if (property_names[i] == "compaction pressure")
else if (property_names[i] == "total pressure")
{
const unsigned int pc_comp_idx = this->introspection().variable("compaction pressure").first_component_index;
const double p_c_bar = input_data.solution_values[q][pc_comp_idx];

computed_quantities[q][output_index] = p_c_scale * p_c_bar;
const unsigned int pc_comp_idx = this->introspection().variable("total pressure").first_component_index;
computed_quantities[q][output_index] = input_data.solution_values[q][pc_comp_idx];
}
else if (property_names[i] == "darcy coefficient")
{
const double K_D = this->get_melt_handler().limited_darcy_coefficient(melt_outputs->permeabilities[q] / melt_outputs->fluid_viscosities[q], p_c_scale > 0);
const double K_D = melt_outputs->permeabilities[q] / melt_outputs->fluid_viscosities[q];
computed_quantities[q][output_index] = K_D;
}
else if (property_names[i] == "darcy coefficient no cutoff")
Expand Down Expand Up @@ -241,8 +234,8 @@ namespace aspect
"'Postprocess/Visualization/Melt material properties/List of properties' contains entries more than once. "
"This is not allowed. Please check your parameter file."));

// Always output compaction pressure
property_names.insert(property_names.begin(),"compaction pressure");
// Always output total pressure
property_names.insert(property_names.begin(),"total pressure");
}
prm.leave_subsection();
}
Expand Down
6 changes: 3 additions & 3 deletions source/simulator/assembly.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -383,7 +383,7 @@ namespace aspect
+ finite_element.base_element(introspection.base_elements.pressure).dofs_per_cell;

if (parameters.include_melt_transport)
stokes_dofs_per_cell += finite_element.base_element(introspection.variable("compaction pressure").base_index).dofs_per_cell;
stokes_dofs_per_cell += finite_element.base_element(introspection.variable("total pressure").base_index).dofs_per_cell;

auto worker = [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
internal::Assembly::Scratch::StokesPreconditioner<dim> &scratch,
Expand Down Expand Up @@ -533,7 +533,7 @@ namespace aspect
std::vector<std::vector<bool>> constant_modes;
dealii::ComponentMask cm_pressure = introspection.component_masks.pressure;
if (parameters.include_melt_transport)
cm_pressure = cm_pressure | introspection.variable("compaction pressure").component_mask;
cm_pressure = cm_pressure | introspection.variable("total pressure").component_mask;
DoFTools::extract_constant_modes (dof_handler,
cm_pressure,
constant_modes);
Expand Down Expand Up @@ -805,7 +805,7 @@ namespace aspect
+ finite_element.base_element(introspection.base_elements.pressure).dofs_per_cell;

if (parameters.include_melt_transport)
stokes_dofs_per_cell += finite_element.base_element(introspection.variable("compaction pressure").base_index).dofs_per_cell;
stokes_dofs_per_cell += finite_element.base_element(introspection.variable("total pressure").base_index).dofs_per_cell;

const bool use_reference_density_profile = (parameters.formulation_mass_conservation == Parameters<dim>::Formulation::MassConservation::reference_density_profile)
|| (parameters.formulation_mass_conservation == Parameters<dim>::Formulation::MassConservation::implicit_reference_density_profile);
Expand Down
30 changes: 15 additions & 15 deletions source/simulator/core.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -983,23 +983,23 @@ namespace aspect
}

coupling[x.velocities[d]][
introspection.variable("compaction pressure").first_component_index] = DoFTools::always;
coupling[introspection.variable("compaction pressure").first_component_index]
[x.velocities[d]]
= DoFTools::always;
coupling[x.velocities[d]]
[introspection.variable("fluid pressure").first_component_index]
= DoFTools::always;
coupling[introspection.variable("fluid pressure").first_component_index]
introspection.variable("total pressure").first_component_index] = DoFTools::always;
coupling[introspection.variable("total pressure").first_component_index]
[x.velocities[d]]
= DoFTools::always;
}

coupling[introspection.variable("fluid pressure").first_component_index]
[introspection.variable("fluid pressure").first_component_index]
= DoFTools::always;
coupling[introspection.variable("compaction pressure").first_component_index]
[introspection.variable("compaction pressure").first_component_index]
coupling[introspection.variable("total pressure").first_component_index]
[introspection.variable("total pressure").first_component_index]
= DoFTools::always;
coupling[introspection.variable("total pressure").first_component_index]
[introspection.variable("fluid pressure").first_component_index]
= DoFTools::always;
coupling[introspection.variable("fluid pressure").first_component_index]
[introspection.variable("total pressure").first_component_index]
= DoFTools::always;
}
else
Expand Down Expand Up @@ -1274,12 +1274,12 @@ namespace aspect
{
coupling[introspection.variable("fluid pressure").first_component_index]
[introspection.variable("fluid pressure").first_component_index] = DoFTools::always;
coupling[introspection.variable("compaction pressure").first_component_index]
[introspection.variable("compaction pressure").first_component_index] = DoFTools::always;
coupling[introspection.variable("compaction pressure").first_component_index]
coupling[introspection.variable("total pressure").first_component_index]
[introspection.variable("total pressure").first_component_index] = DoFTools::always;
coupling[introspection.variable("total pressure").first_component_index]
[introspection.variable("fluid pressure").first_component_index] = DoFTools::always;
coupling[introspection.variable("fluid pressure").first_component_index]
[introspection.variable("compaction pressure").first_component_index] = DoFTools::always;
[introspection.variable("total pressure").first_component_index] = DoFTools::always;
}
else
coupling[x.pressure][x.pressure] = DoFTools::always;
Expand Down Expand Up @@ -1560,7 +1560,7 @@ namespace aspect
{
introspection.index_sets.locally_owned_melt_pressure_dofs = system_index_set & Utilities::extract_locally_active_dofs_with_component(dof_handler,
introspection.variable("fluid pressure").component_mask|
introspection.variable("compaction pressure").component_mask);
introspection.variable("total pressure").component_mask);
introspection.index_sets.locally_owned_fluid_pressure_dofs = system_index_set & Utilities::extract_locally_active_dofs_with_component(dof_handler,
introspection.variable("fluid pressure").component_mask);
}
Expand Down
Loading
Loading