Commit cfa1e637 authored by Udo Ziegler's avatar Udo Ziegler
Browse files

version 4.1 update

parent c3b4600d
/** @file configUser.c
*
* @brief User-defined initial conditions
*
* Copyright (C) 2020 Udo Ziegler
*
* Category: Numerics.User
*
* Version: 4.1
*
* Public functions: configUser()
*/
#include <nirvana.h>
/*############################################################################*/
void configUser(GRD **gm)
/**
* Definition of ICs.
*
* @param gm master mesh pointer
*
* @author Udo Ziegler, AIP
*
* @note Testproblem: MHD/problem29: Linear wave propagation
* (cf. Matsumoto+ 19, ApJ 874, 37)
*/
{
GRD *g;
double g1,gammaa,gammaa2,pert,vx0,a,va,rho0,p0,b0,bx0,by0,lambda,vx1,vy1,
by1,rho1,p1,psin;
int ix,iy,iz;
_C.permeability=1.;
vx0=_C.param[0];
a=_C.param[1];
va=_C.param[2];
g1=_C.gamma-1.;
gammaa=BORIS_CORR ? 1./sqrt(1.+va*va/BORIS_CRED/BORIS_CRED) : 1.;
gammaa2=gammaa*gammaa;
pert=1.e-4;
rho0=1.;
p0=a*a*rho0/_C.gamma;
b0=va*sqrt(rho0);
switch(_C.flag[0])
{
case 1: /* Alfven wave */
bx0=b0;
by0=0.;
lambda=0.5*((1+gammaa2)*vx0+sqrt((1-gammaa2)*(1-gammaa2)*vx0*vx0
+4.*gammaa2*va*va));
vx1=0.;
vy1=lambda-vx0;
by1=-b0;
break;
case 2: /* sound wave */
bx0=b0;
by0=0.;
lambda=0.5*((_C.gamma*(1-gammaa2)+2.*gammaa2)*vx0+sqrt((1-gammaa2)*(1-gammaa2)
*vx0*vx0*_C.gamma*_C.gamma+4.*(_C.gamma-1.)*gammaa2*(1.-gammaa2)
*vx0*vx0+4.*gammaa2*a*a));
rho1=rho0;
vx1=lambda-vx0;
p1=((_C.gamma-1.)*vx0*(vx0-lambda)*(vx0-lambda)+(vx0-lambda)*a*a)
/(_C.gamma*vx0-lambda);
vy1=by1=0.;
break;
case 3: /* magnetosonic wave */
bx0=0.;
by0=b0;
lambda=sqrt(va*va+a*a);
rho1=rho0;
vx1=gammaa*lambda;
vy1=0.;
by1=b0;
p1=_C.gamma*p0;
break;
}
for(g=gm[0]; (g); g=g->next)
{
for(iz=0; iz<=g->nz; iz++)
for(iy=0; iy<=g->ny; iy++)
for(ix=0; ix<=g->nx+1; ix++)
g->bx[iz][iy][ix]=bx0;
for(iz=0; iz<=g->nz; iz++)
for(iy=0; iy<=g->ny+1; iy++)
for(ix=0; ix<=g->nx; ix++)
g->by[iz][iy][ix]=by0+by1*pert*sin(2.*PI*g->xc[ix]);
for(iz=0; iz<=g->nz+_C.idz; iz++)
for(iy=0; iy<=g->ny; iy++)
for(ix=0; ix<=g->nx; ix++)
g->bz[iz][iy][ix]=0.;
for(iz=0; iz<=g->nz; iz++)
for(iy=0; iy<=g->ny; iy++)
for(ix=0; ix<=g->nx; ix++)
{
psin=pert*sin(2.*PI*g->xc[ix]);
g->rho[iz][iy][ix]=rho0+rho1*psin;
g->mx[iz][iy][ix]=(vx0+vx1*psin)*g->rho[iz][iy][ix];
g->my[iz][iy][ix]=vy1*psin*g->rho[iz][iy][ix];
g->mz[iz][iy][ix]=0.;
g->e[iz][iy][ix]=(p0+p1*psin)/g1+0.5/g->rho[iz][iy][ix]
*SPQR(g->mx[iz][iy][ix],g->my[iz][iy][ix],0.)
+SPQR(g->bx[iz][iy][ix],g->by[iz][iy][ix],0.)/2.;
}
}
if(MASTER) fprintf(_FLOG,"..INFO(configUser): MHD/PROBLEM29 INITIALIZED.\n");
return;
}
/** @file nirvana.par
*
* @brief NIRVANA parameter file
*
* Specification of basic code parameter for a problem.
*
* Copyright (C) 2020 Udo Ziegler
*
* Category: Infrastructure.Kernel
*
* Version: 4
*
* @author Udo Ziegler, AIP
*
* @note Testproblem: MHD/problem29: Linear wave propagation
*/
>SIMULATION I/O ----------------------------------------------------------------------------------
NEW ./NIRLAST >mode:{NEW,RES,MOD,ANA},fname
100020 1.00e+00 >mod_max,time_max
0010 00500 91005 3.60e+03 >freq_{log,nir,ana},freq_walltime
>GEOMETRY ----------------------------------------------------------------------------------------
CART 0.00e+00 0.00e+00 0.00e+00 >geometry:{CART,CYL,SPH},omega[0-2]
>DOMAIN SETTINGS ---------------------------------------------------------------------------------
-0.50000000e+00 0.50000000e+00 0128 >lo[0],up[0],dim[0]
0.00000000e+00 1.00000000e+00 0004 >lo[1],up[1],dim[1]
0.00000000e+00 1.00000000e+00 0000 >lo[2],up[2],dim[2]
BLOCK 64 1 0 >partitioning_type:{SFC,BLOCK # # #}
>BOUNDARY CONDITIONS -----------------------------------------------------------------------------
PPPPPP >bc[0-5]:{I,O,A,M,R,C,P,D,F,U}
>MESH REFINEMENT ---------------------------------------------------------------------------------
0 0 >imr,amr
2.0e-01 4.0e-01 3.0e-01 -1.0e-00 2.0e-01 >amr_eps[0-4] / refinement thresholds
0.0e-00 1.0e-01 0.0e+00 1.0e-12 1.0e-02 >rhoref,vref,eref,Bref,Cref
0.6e+00 0.2e+00 >amr_{d1,exp}
2.0e-01 0.0e+00 -0.2e+00 >amr_{Jeans,dJeans,Field}
>USER-SPECIFIC PARAMETERS ------------------------------------------------------------------------
4.70000e+00 1.30000e+00 5.00000e+00 >param[0-2] / vx0,a,va
0.00000e+00 0.00000e+00 0.00000e+00 >param[3-5] /
0.00000e+00 0.00000e+00 0.00000e+00 >param[6-8] /
0.00000e+00 0.00000e+00 0.00000e+00 >param[9-11] /
0.00000e+00 0.00000e+00 0.00000e+00 >param[12-14] /
3 0 0 0 0 0 0 0 >flag[0-7] / 1=Alfven, 2=sound, 3=magnetosonic
>SOLVER SPECIFICATIONS ---------------------------------------------------------------------------
RK3 HLLD CT 4.0e-01 >mhd_solver_{time,flux,ef},mhd_courant
STD 4.0e-01 >viscosity_{solver,courant}
STD 4.0e-01 >diffusion_{solver,courant}
STD 4.0e-01 >conduction_{solver,courant}
STD 4.0e-01 >APdiffusion_{solver,courant}
0.1 >heatloss_max_change
0.1 0.1 >reactions_max_change{X,T}
1. >dt0_reduce:{<=1.}
>PHYSICS SPECIFICATIONS --------------------------------------------------------------------------
Y 7.94e+05 >mhd,permeability_rel
N 0.00e+00 >viscosity{,_coeff}
N 0.00e+00 >diffusion{,_coeff}
N 0.00e+00 0.00e+00 0.00e+00 >conduction{,_coeff,_coeff_perp,_sat}
N 0.00e+00 >APdiffusion{,_coeff}
Y 0.01 >energy,energy_dual_sw:[0,1]
N >force
N >gravity
ADI 1.40e+00 1.00e+00 1.00e+00 >eos,gamma,poly_const,temperature
N >heatloss
0 >tracer
0 >fields
N N N 1.0 0.00 >species,reac,energy_reac,mmw,mean_ionization
>>END INPUT --------------------------------------------------------------------------------------
/** @file nirvanaUser.h
*
* @brief User-controllable NIRVANA header
*
* NIRVANA header file containing user-controllable macros.
*
* Copyright (C) 2020 Udo Ziegler
*
* Category: Infrastructure.User
*
* Version: 4.1
*
* @author Udo Ziegler, AIP
*
* @note DO NOT REMOVE ANY MACRO
*
* Testproblem: MHD/problem29: Linear wave propagation
*/
/* NUMERICS-RELATED PARAMETERS ---------------------------------------------- */
#define SPACE_ORDER 2 /* {2} */
#define LIM(a) VL(a) /* SLOPE LIMITER:{MM,MC,VL} */
#define LIM_MULTI_DIM YES /* {YES,NO} */
#define PSI 0.9 /* PARAMETER IN MULTI_D_LIMITER:[0.5,1] */
#define MAXLEVEL 24 /* MAXIMUM REFINEMENT LEVEL: (<128) */
#define MG_ITMAX 80 /* MULTI-GRID: MAX ITERATIONS */
#define MG_TYPE MULT /* MULTI-GRID: TYPE:{MULT} */
#define MG_TOL 1.e-6 /* MULTI-GRID: ERROR TOLERANCE */
#define RKL_COURANT_EXPL 0.25 /* RKL: EXPL COURANT NUMBER: <0.5 */
#define RKL_MAX_COURANT 1.e3 /* RKL: MAX COURANT NUMBER: <=1000 */
#define RKL_DT_LIM 0.10 /* RKL: DT/DT_DYN LIMIT */
#define HEATLOSS_TOL 1.e-4 /* HEATLOSS SOLVER: REL. ERR TOL */
#define HEATLOSS_ATOL 0.01 /* HEATLOSS SOLVER: ABS. ERR TOL */
#define REACTIONS_TOL 1.e-5 /* REACTIONS SOLVER: REL. ERR TOL */
#define REACTIONS_ATOL_X 0. /* REACTIONS SOLVER: ABS. ERR TOL_X */
#define REACTIONS_ATOL_T 0.01 /* REACTIONS SOLVER: ABS. ERR TOL_T */
#define HLLD_PRESSURE_CORR NO /* HLLD: LOW MACH CORRECTION:{NO,YES} */
#define BORIS_CORR YES /* BORIS CORRECTION:{NO,YES,AUTO} */
#define BORIS_CRED 10.0 /* BORIS: RED. C */
#define BORIS_CRED_MIN -1. /* BORIS: MIN RED. C IN UNITS 2|v|+cs */
#define BORIS_AUTO_CA_MAX 3. /* BORIS: MAX CA IN UNITS 2|v|+cs */
/* PHYSICS-RELATED PARAMETERS ----------------------------------------------- */
#define CENTRIFUGAL_FORCE YES /* CENTRIFUGAL FORCE TERM:{YES,NO} */
#define COND_FORCE_ISO NO /* ISOTROPY IN CONDUCTION:{YES,NO} */
#define B_FIELD_SPLITTING NO /* MAGNETIC FIELD SPLITTING:{YES,NO} */
/* NCCM-RELATED PARAMETERS -------------------------------------------------- */
#define CR_IH 1.75e-16 /* H CR IONISATION RATE STD: 2.5e-17
Indriolo & McCall: 1.75e-16 */
#define TDUST 10. /* DUST TEMPERATURE */
#define ADUST 1. /* DUST ABUNDANCE REL. TO MILKY WAY */
#define ISRF_G 0.03 /* ISRF HABING STRENGTH (STD=1.13) */
#define DBNCCM YES /* NCCM DATA BASE AVAILABILITY FLAG */
#define CHECK_REACTIONS YES /* CHECK REACTIONS:{YES,NO} */
/* ANALYTIC EOS-RELATED MACROS ---------------------------------------------- */
#define PUSR(rho,eth) (0.*(rho)+0.*(eth))
#define CS2USR(rho,eth) (0.*(rho)+0.*(eth))
#define TUSR(rho,eth) (0.*(rho)+0.*(eth))
#define ETUSR(rho,eth) (0.*(rho)+0.*(eth))
/** @file analysisUser.c
*
* @brief User-defined data analysis and runtime interaction
*
* Copyright (C) 2020 Udo Ziegler
*
* Category: Numerics.User
*
* Version: 4
*
* Public functions: analysisUser()
*/
#include <nirvana.h>
/*############################################################################*/
void analysisUser(GRD **gm)
/**
* Template for user-defined data analysis and runtime interaction tasks.
*
* @param gm master mesh pointer
*
* @author Udo Ziegler, AIP
*
* @note Testproblem: MHD/problem30: Magnetosphere-wind interaction
*/
{
FILE *fpt;
GRD *g;
double rhow,Ti,ei,Ri,Ri2,bx,by,bz,pih,dy;
int il,ix,iy,iz;
char fname[BUFS];
/* MY IMPLEMENTATION */
/* SET INNER BOUNDARY */
rhow=_C.param[5]*_C.mean_molecular_weight*AMU;
Ti=3.5e+04;
ei=_C.param[5]*KB*Ti/(_C.gamma-1.);
Ri=_C.param[1]*RE;
pih=0.5/_C.permeability;
Ri2=Ri*Ri;
for(il=0; il<=_C.level_all; il++)
for(g=gm[il]; (g); g=g->next)
{
for(iz=0; iz<=g->nz; iz++)
for(iy=0; iy<=g->ny; iy++)
for(ix=0; ix<=g->nx; ix++)
if(SPQR(g->xc[ix],g->yc[iy],g->zc[iz])<Ri2)
{
g->rho[iz][iy][ix]=rhow;
g->mx[iz][iy][ix]=g->my[iz][iy][ix]=g->mz[iz][iy][ix]=0.;
}
for(iz=g->izs; iz<=g->ize; iz++)
for(iy=g->iys; iy<=g->iye; iy++)
for(ix=g->ixs; ix<=g->ixe; ix++)
if(SPQR(g->xc[ix],g->yc[iy],g->zc[iz])<Ri2)
{
if(SPQR(g->xc[ix+1],g->yc[iy],g->zc[iz])>Ri2)
g->mx[iz][iy][ix]=-g->mx[iz][iy][ix+1]*rhow/g->rho[iz][iy][ix+1];
if(SPQR(g->xc[ix-1],g->yc[iy],g->zc[iz])>Ri2)
g->mx[iz][iy][ix]=-g->mx[iz][iy][ix-1]*rhow/g->rho[iz][iy][ix-1];
if(SPQR(g->xc[ix],g->yc[iy+1],g->zc[iz])>Ri2)
g->my[iz][iy][ix]=-g->my[iz][iy+1][ix]*rhow/g->rho[iz][iy+1][ix];
if(SPQR(g->xc[ix],g->yc[iy-1],g->zc[iz])>Ri2)
g->my[iz][iy][ix]=-g->my[iz][iy-1][ix]*rhow/g->rho[iz][iy-1][ix];
if(SPQR(g->xc[ix],g->yc[iy],g->zc[iz+1])>Ri2)
g->mz[iz][iy][ix]=-g->mz[iz+1][iy][ix]*rhow/g->rho[iz+1][iy][ix];
if(SPQR(g->xc[ix],g->yc[iy],g->zc[iz-1])>Ri2)
g->mz[iz][iy][ix]=-g->mz[iz-1][iy][ix]*rhow/g->rho[iz-1][iy][ix];
}
for(iz=g->izs; iz<=g->ize; iz++)
for(iy=g->iys; iy<=g->iye; iy++)
for(ix=g->ixs; ix<=g->ixe; ix++)
if(SPQR(g->xc[ix],g->yc[iy],g->zc[iz])<Ri2)
{
dy=(g->yc[iy]-g->y[iy])*g->dyi;
bx=g->bx[iz][iy][ix]+(g->bx[iz][iy][ix+1]-g->bx[iz][iy][ix])
*(g->xc[ix]-g->x[ix])*g->dxi;
by=g->by[iz][iy][ix]+(g->by[iz][iy+1][ix]-g->by[iz][iy][ix])*dy;
bz=0.5*(g->bz[iz][iy][ix]+g->bz[iz+1][iy][ix]);
g->e[iz][iy][ix]=ei+pih*SPQR(bx,by,bz)
+0.5*SPQR(g->mx[iz][iy][ix],g->my[iz][iy][ix],g->mz[iz][iy][ix])/rhow;
}
}
return;
}
/** @file bcbUser.c
*
* @brief User-defined boundary conditions for the magnetic field
*
* Copyright (C) 2020 Udo Ziegler
*
* Category: Numerics.User
*
* Version: 4
*
* Public functions: bcbUser()
*/
#include <nirvana.h>
/*############################################################################*/
void bcbUser(GRD *g, flag_t boundary)
/**
* Template for user-defined boundary conditions for the magnetic field.
*
* @param g grid pointer
* @param boundary domain boundary:{X_LOWER,X_UPPER,Y_LOWER,Y_UPPER,
* Z_LOWER,Z_UPPER}
*
* @author Udo Ziegler, AIP
*
* @note Testproblem: MHD/problem30: Magnetosphere-wind interaction
*/
{
double bzi;
int ix,iy,iz;
bzi=_C.param[3];
switch(boundary)
{
case X_LOWER: /* LOWER X-FACE --------------------------------------------- */
for(ix=0; ix<g->ixs; ix++)
for(iz=0; iz<=g->nz; iz++)
for(iy=0; iy<=g->ny; iy++)
{
}
break;
case X_UPPER: /* UPPER X-FACE --------------------------------------------- */
for(ix=g->ixe+1; ix<=g->nx; ix++)
for(iz=0; iz<=g->nz; iz++)
for(iy=0; iy<=g->ny; iy++)
{
g->bx[iz][iy][ix+1]=0.;
g->by[iz][iy][ix]=0.;
g->bz[iz][iy][ix]=bzi;
}
break;
case Y_LOWER: /* LOWER Y-FACE --------------------------------------------- */
for(iy=0; iy<g->iys; iy++)
for(iz=0; iz<=g->nz; iz++)
for(ix=0; ix<=g->nx; ix++)
{
}
break;
case Y_UPPER: /* UPPER Y-FACE --------------------------------------------- */
for(iy=g->iye+1; iy<=g->ny; iy++)
for(iz=0; iz<=g->nz; iz++)
for(ix=0; ix<=g->nx; ix++)
{
}
break;
case Z_LOWER: /* LOWER Z-FACE --------------------------------------------- */
for(iz=0; iz<g->izs; iz++)
for(iy=0; iy<=g->ny; iy++)
for(ix=0; ix<=g->nx; ix++)
{
}
break;
case Z_UPPER: /* UPPER Z-FACE --------------------------------------------- */
for(iz=g->ize+1; iz<=g->nz; iz++)
for(iy=0; iy<=g->ny; iy++)
for(ix=0; ix<=g->nx; ix++)
{
}
break;
}
return;
}
/** @file bceUser.c
*
* @brief User-defined boundary conditions for the total energy density.
*
* Copyright (C) 2020 Udo Ziegler
*
* Category: Numerics.User
*
* Version: 4
*
* Public functions: bceUser()
*/
#include <nirvana.h>
/*############################################################################*/
void bceUser(GRD *g, flag_t boundary)
/**
* Template for user-defined boundary conditions for the total energy
* density.
*
* @param g grid pointer
* @param boundary domain boundary:{X_LOWER,X_UPPER,Y_LOWER,Y_UPPER,
* Z_LOWER,Z_UPPER}
*
* @author Udo Ziegler, AIP
*
* @note Testproblem: MHD/problem30: Magnetosphere-wind interaction
*/
{
double ew;
int ix,iy,iz;
ew=_C.param[5]*KB*_C.param[6]/(_C.gamma-1.)
+0.5*_C.mean_molecular_weight*AMU*_C.param[5]*_C.param[4]*_C.param[4]
+0.5*_C.param[3]*_C.param[3]/_C.permeability;
switch(boundary)
{
case X_LOWER: /* LOWER X-FACE --------------------------------------------- */
for(ix=0; ix<g->ixs; ix++)
for(iz=0; iz<=g->nz; iz++)
for(iy=0; iy<=g->ny; iy++)
{
}
break;
case X_UPPER: /* UPPER X-FACE --------------------------------------------- */
for(ix=g->ixe+1; ix<=g->nx; ix++)
for(iz=0; iz<=g->nz; iz++)
for(iy=0; iy<=g->ny; iy++)
{
g->e[iz][iy][ix]=ew;
}
break;
case Y_LOWER: /* LOWER Y-FACE --------------------------------------------- */
for(iy=0; iy<g->iys; iy++)
for(iz=0; iz<=g->nz; iz++)
for(ix=0; ix<=g->nx; ix++)
{
}
break;
case Y_UPPER: /* UPPER Y-FACE --------------------------------------------- */
for(iy=g->iye+1; iy<=g->ny; iy++)
for(iz=0; iz<=g->nz; iz++)
for(ix=0; ix<=g->nx; ix++)
{
}
break;
case Z_LOWER: /* LOWER Z-FACE --------------------------------------------- */
for(iz=0; iz<g->izs; iz++)
for(iy=0; iy<=g->ny; iy++)
for(ix=0; ix<=g->nx; ix++)
{
}
break;
case Z_UPPER: /* UPPER Z-FACE --------------------------------------------- */
for(iz=g->ize+1; iz<=g->nz; iz++)
for(iy=0; iy<=g->ny; iy++)
for(ix=0; ix<=g->nx; ix++)
{
}
break;
}
return;
}
/** @file bcetUser.c
*
* @brief User-defined boundary conditions for the thermal energy density.
*
* Copyright (C) 2020 Udo Ziegler
*
* Category: Numerics.User
*
* Version: 4
*
* Public functions: bcetUser()
*/
#include <nirvana.h>
/*############################################################################*/
void bcetUser(GRD *g, flag_t boundary)
/**
* Template for user-defined boundary conditions for the thermal energy
* density.
*
* @param g grid pointer
* @param boundary domain boundary-{X_LOWER,X_UPPER,Y_LOWER,Y_UPPER,
* Z_LOWER,Z_UPPER}
*
* @author Udo Ziegler, AIP
*
* @note Testproblem: MHD/problem30: Magnetosphere-wind interaction
*/
{
double ew;