Page MenuHomec4science
Files F83529437

art_visc.c
No OneTemporary
Actions

File Metadata

Created
Tue, Sep 17, 15:36

art_visc.c
View Options

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_matrix.h>
#include "allvars.h"
#include "proto.h"
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_ROSS) || defined(ART_VISCO_CD)
#ifdef ART_VISCO_CD
int sgsl;
gsl_matrix * Dgsl;
gsl_matrix * Tgsl;
gsl_matrix * Tigsl;
gsl_matrix * Vgsl;
gsl_matrix * Vtgsl;
gsl_matrix * Idgsl;
gsl_permutation * perm;
double DiVelAccurate;
double dt_DiVel;
double shock_indicator_CD, shock_limiter_CD;
double h_i2, MaxSignalVel2;
double alpha_CD_max;
double tr_S;
double Xi_num, Xi_den;
void art_visc_allocate()
{
Dgsl = gsl_matrix_alloc (3, 3);
Tgsl = gsl_matrix_alloc (3, 3);
Tigsl = gsl_matrix_alloc (3, 3);
perm = gsl_permutation_alloc (3);
Vgsl = gsl_matrix_alloc (3, 3);
Vtgsl = gsl_matrix_alloc (3, 3);
Idgsl = gsl_matrix_alloc (3, 3);
}
void art_visc_free()
{
gsl_matrix_free (Dgsl);
gsl_matrix_free (Tgsl);
gsl_matrix_free (Tigsl);
gsl_matrix_free (Vgsl);
gsl_matrix_free (Vtgsl);
gsl_matrix_free (Idgsl);
gsl_permutation_free (perm);
}
void compute_art_visc(int i)
{
int ii,jj;
/* FILL GSL MATRIX*/
for (ii = 0; ii < 3; ii++)
for (jj = 0; jj < 3; jj++)
{
gsl_matrix_set (Dgsl, ii, jj, SphP[i].DmatCD[ii][jj]);
gsl_matrix_set (Tgsl, ii, jj, SphP[i].TmatCD[ii][jj]);
}
/* COMPUTE INV(TmatCD) */
// Make LU decomposition of matrix Tgsl
gsl_linalg_LU_decomp (Tgsl, perm, &sgsl);
// Invert the matrix Tgsl
gsl_linalg_LU_invert (Tgsl, perm, Tigsl);
/* DiVelAccurate = COMPUTE Tr(D * INV(TmatCD)) */
gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
1.0, Dgsl, Tigsl,
0.0, Vgsl);
/* new div(velpred) */
DiVelAccurate = gsl_matrix_get(Vgsl, 0, 0) + gsl_matrix_get(Vgsl, 1, 1) + gsl_matrix_get(Vgsl, 2, 2);
/* d/dt(div(velpred)) */
dt_DiVel = (P[i].Ti_endstep - P[i].Ti_begstep) * All.Timebase_interval;
if(All.ComovingIntegrationOn)
{
printf("ART_VISCO_CD doesnot work in cosmo! \n");
endrun(768543);
}
/* !!!!!! for the initial conditions !!!!! */
if (dt_DiVel>0.)
SphP[i].DiVelTemp = (DiVelAccurate - SphP[i].DiVelAccurate) / dt_DiVel;
else
SphP[i].DiVelTemp = 0.;
/* compute chock limiter R_i, Xi_i*/
SphP[i].R_CD = SphP[i].R_CD / SphP[i].Density;
for (ii = 0; ii < 3; ii++)
for (jj = 0; jj < 3; jj++)
{
gsl_matrix_set( Vtgsl, jj, ii, gsl_matrix_get( Vgsl, ii, jj ) );
}
gsl_matrix_set_identity (Idgsl);
gsl_matrix_add (Vgsl, Vtgsl);
gsl_matrix_scale (Vgsl, .5);
gsl_matrix_scale (Idgsl, -SphP[i].DiVelAccurate/3.);
/* tensor S */
gsl_matrix_add (Vgsl, Idgsl );
for (ii = 0; ii < 3; ii++)
for (jj = 0; jj < 3; jj++)
{
gsl_matrix_set( Vtgsl, jj, ii, gsl_matrix_get( Vgsl, ii, jj ) );
}
gsl_matrix_mul_elements(Vgsl, Vtgsl);
tr_S = gsl_matrix_get(Vgsl, 0, 0) + gsl_matrix_get(Vgsl, 1, 1) + gsl_matrix_get(Vgsl, 2, 2);
Xi_num = pow(2. * pow(1. - SphP[i].R_CD, 4.) * SphP[i].DiVelAccurate, 2);
Xi_den = Xi_num + tr_S ;
/* update */
SphP[i].DiVelAccurate = DiVelAccurate;
// printf("%4.2f %4.2f \n",DiVelAccurate, SphP[i].DivVel );
alpha_CD_max = All.ArtBulkViscConstMax;
/*shock_limiter_CD = 1. ;*/
shock_limiter_CD = Xi_num / Xi_den;
shock_indicator_CD = shock_limiter_CD * dmax(-SphP[i].DiVelTemp,0.);
h_i2 = SphP[i].Hsml;
h_i2 = h_i2 * h_i2;
MaxSignalVel2 = SphP[i].MaxSignalVel;
MaxSignalVel2 = MaxSignalVel2 * MaxSignalVel2;
SphP[i].ArtBulkViscConstOld = alpha_CD_max * h_i2 * shock_indicator_CD / (MaxSignalVel2 + h_i2 * shock_indicator_CD);
// printf("%4.2f \n",SphP[i].alpha_loc);
// exit(1);
}
#endif
void move_art_visc(int i,double dt_drift)
{
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_ROSS)
double timedecay;
double source_Velbased;
double p_over_rho2_i, soundspeed_i;
double dalpha_over_dt;
#endif
#ifdef ART_VISCO_CD
double timedecay, dalpha_over_dt;
#endif
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_ROSS)
source_Velbased = dmax(-SphP[i].DivVel,0.) ;
p_over_rho2_i = SphP[i].Pressure / (SphP[i].Density * SphP[i].Density);
soundspeed_i = sqrt(GAMMA * p_over_rho2_i * SphP[i].Density);
timedecay = SphP[i].Hsml/(2. * All.ArtBulkViscConstL * soundspeed_i);
#ifdef ART_VISCO_MM
dalpha_over_dt = (All.ArtBulkViscConstMin - SphP[i].alpha) / timedecay + source_Velbased;
#endif
#ifdef ART_VISCO_ROSS
dalpha_over_dt = (All.ArtBulkViscConstMin - SphP[i].alpha) / timedecay + (All.ArtBulkViscConstMax - SphP[i].alpha) * source_Velbased;
#endif
SphP[i].ArtBulkViscConst = SphP[i].ArtBulkViscConst + dalpha_over_dt * dt_drift;
#endif
#ifdef ART_VISCO_CD
timedecay = SphP[i].Hsml/(2. * All.ArtBulkViscConstL * SphP[i].MaxSignalVel);
if (SphP[i].ArtBulkViscConst < SphP[i].ArtBulkViscConstOld)
SphP[i].ArtBulkViscConst = SphP[i].ArtBulkViscConstOld;
else
SphP[i].ArtBulkViscConst = SphP[i].ArtBulkViscConstOld + (SphP[i].ArtBulkViscConst - SphP[i].ArtBulkViscConstOld) * exp(-dt_drift/timedecay);
#endif
}
#endif

Event Timeline

c4science · Help