Implement gravitational separation flux

[djbower]

The gravitational separation flux quantifies the energy transport due to melt-solid separation.

In the original C-version of Spider the flux is coded in terms of entropy (code given below), but in Aragog this should be converted to temperature, with additional care taken to ensure correct signs for radial derivatives and gravity (which could be different between the C version and Aragog). As part of this effort, permeability laws could also be encapsulated into their own class structure. A temperature form of the gravitational flux can be found in the original work of Abe (1993), as well as subsequent papers (1995, 1997).

There is a placeholder property to implement this energy flux in solver.py (def gravitational_separation_flux):

    @property
    def gravitational_separation_flux(self) -> np.ndarray:
        """Gravitational separation"""
        raise NotImplementedError

Original C code snippet, and other functions (such as the permeability laws) can be found in the C code:

static PetscErrorCode append_Jgrav( Ctx *E )
{
    /* gravitational separation heat flux at basic nodes */

    PetscErrorCode ierr;
    PetscScalar    Jgrav;
    PetscInt       i,ilo_b,ihi_b,w_b;
    Solution       *S = &E->solution;

    PetscFunctionBeginUser;

    /* loop over all basic internal nodes */
    ierr = DMDAGetCorners(E->da_b,&ilo_b,0,0,&w_b,0,0);CHKERRQ(ierr);
    ihi_b = ilo_b + w_b;

    for(i=ilo_b; i<ihi_b; ++i){
        Jgrav = GetGravitationalHeatFlux( E, &i );
        ierr = VecSetValue(S->Jgrav,i,Jgrav,INSERT_VALUES);CHKERRQ(ierr);
    }

    ierr = VecAssemblyBegin(S->Jgrav);CHKERRQ(ierr);
    ierr = VecAssemblyEnd(S->Jgrav);CHKERRQ(ierr);

    ierr = VecAXPY( S->Jtot, 1.0, S->Jgrav ); CHKERRQ(ierr);

    PetscFunctionReturn(0);
}

PetscScalar GetGravitationalHeatFlux( Ctx *E, PetscInt * ind_ptr )
{
    PetscErrorCode  ierr;
    PetscScalar     porosity,cond1,cond2,F,dv,pres,rho,Sliq,Ssol,phi,Jgrav,temp,phi_s,pres_s,Sval,gphi,smth;
    PetscInt        should_be_two;
    Solution const  *S = &E->solution;
    Mesh const      *M = &E->mesh;
    Parameters const P = E->parameters;
    EOS              *sub_eos;
    EOSEvalData     eval_liq, eval_sol;
    PetscInt        ilo_b,ihi_b,w_b;

    ierr = DMDAGetCorners(E->da_b,&ilo_b,0,0,&w_b,0,0);CHKERRQ(ierr);
    ihi_b = ilo_b + w_b; // this is one more than last index of basic array

    /* no gravitational separation at boundaries (since no mass transfer possible) */
    if( (!*ind_ptr) || (*ind_ptr == ihi_b-1) ){
        return 0.0;
    }

    /* Jgrav requires two phases */
    ierr = EOSCompositeGetSubEOS(P->eos, &sub_eos, &should_be_two);CHKERRQ(ierr);
    if (should_be_two!=2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Expecting two sub-EOSs");
    const EOS Ep0 = sub_eos[0];
    const EOS Ep1 = sub_eos[1];

    ierr = VecGetValues(M->pressure_b,1,ind_ptr,&pres);CHKERRQ(ierr);
    ierr = VecGetValues(S->rho,1,ind_ptr,&rho);CHKERRQ(ierr);
    ierr = VecGetValues(S->phi,1,ind_ptr,&phi);CHKERRQ(ierr);

    /* smoothing based on whether cell (staggered node below) actually has
       melt.  Otherwise, non-physical cooling results */
    ierr = VecGetValues(M->pressure_s,1,ind_ptr,&pres_s);CHKERRQ(ierr);
    ierr = VecGetValues(S->phi_s,1,ind_ptr,&phi_s);CHKERRQ(ierr);
    ierr = VecGetValues(S->S_s,1,ind_ptr,&Sval);CHKERRQ(ierr);

    ierr = EOSGetPhaseBoundary( Ep0, pres, &Sliq, NULL );CHKERRQ(ierr);
    ierr = EOSEval(Ep0, pres, Sliq, &eval_liq);CHKERRQ(ierr);
    ierr = EOSGetPhaseBoundary( Ep1, pres, &Ssol, NULL );CHKERRQ(ierr);
    ierr = EOSEval(Ep1, pres, Ssol, &eval_sol);CHKERRQ(ierr);

    porosity = (eval_sol.rho - rho) / ( eval_sol.rho - eval_liq.rho );

    switch( P->SEPARATION ){
        case 1:
            /* Abe formulation */
            /* these switches depend on the functions below, and are constructed to
               ensure F is a smooth function of porosity.  They are given in Abe in
               terms of the volume fraction of solid, hence the 1.0 minus in the if
               statement.  They depend on the choice of constants to the flow laws */
            /* See Eq. 44 in Abe (1995).  But below we use permeability directly to
               stay connected to the physics */
            cond1 = 0.7714620383592684;
            cond2 = 0.0769618;

            /* solid_volume < cond1 (Abe) is porosity > 1-cond1 (here) */
            if(porosity > cond1){
                /* Stokes settling factor with grainsize squared */
                F = (2.0/9.0) * PetscPowScalar( P->grain, 2.0);
            }
            else if(porosity < cond2){
                /* permeability includes grainsize squared */
                F = GetPermeabilityBlakeKozenyCarman( P->grain, porosity, 1.0E-3 );
                F /= porosity;
            }
            else{
                /* permeability includes grainsize squared */
                F = GetPermeabilityRumpfGupte( P->grain, porosity, 5.0/7.0 );
                F /= porosity;
            }
            break;
        case 2:
            /* numerically it seems problematic to have no separation and then
               try and turn it on for low melt fraction.  Hence this seems
               to work less well than the Abe smooth approach above */
            /* apply separation below 40% melt volume fraction */
            //if( porosity < 0.4 ){
                /* large permeability from Rudge (2018) */
            F = GetPermeabilityRudge( P->grain, porosity, 1.0/75 );
            F /= porosity;
            //}
            //else{
            //    F = 0.0;
            //}
            {
                PetscScalar matprop_smooth_width;
                ierr = EOSCompositeGetMatpropSmoothWidth(P->eos, &matprop_smooth_width);CHKERRQ(ierr);
                F *= 1.0 - tanh_weight( porosity, 0.3, 1.0E-2 );
            }
            break;
        default:
            SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Unsupported SEPARATION value %d provided",P->SEPARATION);
            break;
    }

    /* relative velocity: velocity_melt - velocity_solid */
    /* works for Stokes and permeability since F is different above */
    /* e.g. Abe (1995) Eq. 39 and 40 */
    dv = ( eval_liq.rho - eval_sol.rho ) * P->gravity * F;
    dv /= PetscPowScalar(10.0, P->eos_phases[0]->log10visc);

    /* mass flux, e.g. Abe (1995) Eq. 8 */
    /* here, clear that Jgrav is zero when phi is single phase (phi=0
       or phi=1) */
    Jgrav = rho * phi * ( 1.0 - phi ) * dv;

    /* temperature of the fusion curve */
    temp = eval_sol.T + 0.5 * (eval_liq.T - eval_sol.T);

    /* energy flux */
    Jgrav *= temp * ( Sliq - Ssol ); // enthalpy

    /* smoothing across phase boundaries for two phase composite */
    /* this smooths based on the value of the melt fraction in the cell below
       the interface, but this is only appropriate for bottom-up crystallisation */
    if( P->JGRAV_BOTTOM_UP){
        ierr = EOSCompositeGetTwoPhasePhaseFractionNoTruncation(P->eos, pres_s, Sval, &gphi);CHKERRQ(ierr);
        {
          PetscScalar matprop_smooth_width;

          ierr = EOSCompositeGetMatpropSmoothWidth(P->eos, &matprop_smooth_width);CHKERRQ(ierr);
          smth = get_smoothing(matprop_smooth_width, gphi );
        }

        Jgrav *= smth;
    }

    return Jgrav;
}

[djbower]

Some permeability laws:

static PetscScalar GetPermeabilityBlakeKozenyCarman( PetscScalar grainsize, PetscScalar porosity, PetscScalar constant )
{
    /* Abe (1995) Eq. 42a.  Also see McKenzie (1984) */

    PetscScalar kp;

    kp = PetscPowScalar( grainsize, 2.0 ) * PetscPowScalar( porosity, 3.0 ) / PetscPowScalar( 1.0-porosity, 2.0 );
    kp *= constant;

    return kp;

}

static PetscScalar GetPermeabilityRumpfGupte( PetscScalar grainsize, PetscScalar porosity, PetscScalar constant )
{
    /* Abe (1995) Eq. 42b.  Also see McKenzie (1984) */

    PetscScalar kp;

    kp = PetscPowScalar( grainsize, 2.0 ) * PetscPowScalar( porosity, 5.5 );
    kp *= constant;

    return kp;

}

static PetscScalar GetPermeabilityRudge( PetscScalar grainsize, PetscScalar porosity, PetscScalar constant )
{
    /* Rudge (2018) Eq. 7.4, large permeability for 0.1 < porosity < 0.3 */

    PetscScalar kp;

    kp = PetscPowScalar( grainsize, 2.0 ) * PetscPowScalar( porosity, 3.0 );
    kp *= constant;

    return kp;

}