Homec4science
Diffusion Gear (master) Browse

Gear/src/art_visc.c
b87f34a0ce7emaster

Last Change
Raw File

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_RO) || 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);
double det_gsl;
gsl_linalg_LU_decomp (Tgsl, perm, &sgsl);
det_gsl = gsl_linalg_LU_det (Tgsl, &sgsl);
if (abs(det_gsl)>0.000001)
{
// Invert the matrix Tgsl
gsl_linalg_LU_invert (Tgsl, perm, Tigsl);
}
else
{
for (ii = 0; ii < 3; ii++)
for (jj = 0; jj < 3; jj++)
{
gsl_matrix_set (Tigsl, ii, jj, 0.);
}
}
/* 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;
if (abs(Xi_den)>0.000001)
shock_limiter_CD = Xi_num / Xi_den;
else
shock_limiter_CD = 0.;
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].MaxSignalVelCD;
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_RO)
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_RO)
source_Velbased = dmax(-SphP[i].DivVel,0.) ;
#ifdef DENSITY_INDEPENDENT_SPH
p_over_rho2_i = SphP[i].Pressure / (SphP[i].EgyWtDensity * SphP[i].EgyWtDensity);
soundspeed_i = sqrt(GAMMA * p_over_rho2_i * SphP[i].EgyWtDensity);
#else
p_over_rho2_i = SphP[i].Pressure / (SphP[i].Density * SphP[i].Density);
soundspeed_i = sqrt(GAMMA * p_over_rho2_i * SphP[i].Density);
#endif
timedecay = SphP[i].Hsml/(2. * All.ArtBulkViscConstL * soundspeed_i);
#ifdef ART_VISCO_MM
dalpha_over_dt = (All.ArtBulkViscConstMin - SphP[i].ArtBulkViscConst) / timedecay + source_Velbased;
#endif
#ifdef ART_VISCO_RO
dalpha_over_dt = (All.ArtBulkViscConstMin - SphP[i].ArtBulkViscConst) / timedecay + dmax(-SphP[i].DivVel*(All.ArtBulkViscConstMax - SphP[i].ArtBulkViscConst) ,0.) ;
#endif
SphP[i].ArtBulkViscConst = SphP[i].ArtBulkViscConst + dalpha_over_dt * dt_drift;
SphP[i].ArtBulkViscConst = dmax(SphP[i].ArtBulkViscConst,All.ArtBulkViscConstMin);
SphP[i].ArtBulkViscConst = dmin(SphP[i].ArtBulkViscConst,All.ArtBulkViscConstMax);
#endif
#ifdef ART_VISCO_CD
timedecay = SphP[i].Hsml/(2. * All.ArtBulkViscConstL * SphP[i].MaxSignalVelCD);
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
c4science · Help