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

directory /nirvana/testproblems/MHD/problem24B added

parent cafb702f
/** @file analysisUser.c
*
* @brief User-defined analysis and 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 analysis and interaction purposes.
*
* @param gm master mesh pointer
*
* @author Udo Ziegler, AIP
*
* @note Testproblem: MHD/problem24B: Kelvin-Helmholtz instability (Low Mach)
*/
{
FILE *fpt;
GRD *g;
double *vy,*vk,ak,bk,ck;
int ix,iy,iz,ik;
static int nrows=1;
char fname[99];
/* COMPUTATION OF FOURIER AMPLITUDES */
vy=Array(_C.dim[0]-1);
vk=Array(_C.dim[0]-1);
for(g=gm[0]; (g); g=g->next)
{
iy=(int)(-g->y[_C.ngc]/g->dy+0.5)+_C.ngc-1;
iz=0;
for(ix=g->ixs; ix<=g->ixe; ix++)
vy[ix-g->ixs]=0.5*(g->my[iz][iy][ix]/g->rho[iz][iy][ix]
+g->my[iz][iy+1][ix]/g->rho[iz][iy+1][ix]);
}
for(ik=0; ik<_C.dim[0]; ik++)
{
ak=0.;
bk=0.;
ck=ik*2.*PI/_C.dim[0];
for(ix=0; ix<_C.dim[0]; ix++)
{
ak+=vy[ix]*cos(ix*ck);
bk+=vy[ix]*sin(ix*ck);
}
vk[ik]=sqrt(ak*ak+bk*bk)/_C.dim[0];
}
free(vy);
if(MASTER)
{
sprintf(fname,"%s/MHD_Problem24B.csv",_IODIR);
if(nrows==1)
{
fpt=fopen(fname,"a");
fprintf(fpt,"ncols:%3d nrows:%9d\n",_C.dim[0]+1,nrows);
fprintf(fpt," time Fourier_amplitudes\n");
}
else
{
fpt=fopen(fname,"r+");
rewind(fpt);
fprintf(fpt,"ncols:%3d nrows:%9d\n",_C.dim[0]+1,nrows);
fclose(fpt);
fpt=fopen(fname,"a");
}
nrows++;
fprintf(fpt,"%12.4e ",_C.time);
for(ik=0; ik<_C.dim[0]; ik++)
fprintf(fpt,"%12.4e ",vk[ik]);
fprintf(fpt,"\n");
fclose(fpt);
}
free(vk);
return;
}
/** @file configUser.c
*
* @brief User-defined initial conditions
*
* Copyright (C) 2020 Udo Ziegler
*
* Category: Numerics.User
*
* Version: 4
*
* 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/problem24B: Kelvin-Helmholtz instability
* (Minoshima+ 2020, ApJS 248:12)
*/
{
GRD *g;
double x,y,lambda,rho0,p0,V0,e0,B0,theta,pi8;
int ix,iy,iz;
if(_C.mf)
{
_C.permeability=1.;
pi8=0.125/_C.permeability;
}
lambda=1.;
rho0=1.;
p0=50000.;
V0=1.;
e0=p0/(_C.gamma-1.);
B0=1.;
theta=71.565/360.*2.*PI;
for(g=gm[0]; (g); g=g->next)
{
if(_C.mf) /* MAGNETIC FIELD */
{
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]=B0*cos(theta);
for(iy=0; iy<=g->ny+1; iy++)
for(ix=0; ix<=g->nx; ix++)
g->by[iz][iy][ix]=0.;
}
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]=B0*sin(theta);
}
/* HD VARIABLES */
for(iz=0; iz<=g->nz; iz++)
for(iy=0; iy<=g->ny; iy++)
for(ix=0; ix<=g->nx; ix++)
{
x=g->xc[ix];
y=g->yc[iy];
g->rho[iz][iy][ix]=rho0;
g->mx[iz][iy][ix]=0.5*rho0*V0*tanh(y/lambda);
//g->my[iz][iy][ix]=fabs(y)<g->dy ? rho0*0.01*V0*sin(2.*PI*x/20./lambda) : 0.;
g->my[iz][iy][ix]=fabs(y)<lambda ? rho0*0.01*V0*sin(2.*PI*x/20./lambda) : 0.;
g->mz[iz][iy][ix]=0.;
g->e[iz][iy][ix]=e0+0.5*SPQR(g->mx[iz][iy][ix],g->my[iz][iy][ix],
g->mz[iz][iy][ix])/rho0;
if(_C.mf) g->e[iz][iy][ix]+=pi8*((g->bx[iz][iy][ix]+g->bx[iz][iy][ix+1])
*(g->bx[iz][iy][ix]+g->bx[iz][iy][ix+1])
+(g->by[iz][iy][ix]+g->by[iz][iy+1][ix])
*(g->by[iz][iy][ix]+g->by[iz][iy+1][ix])
+(g->bz[iz][iy][ix]+g->bz[iz+_C.idz][iy][ix])
*(g->bz[iz][iy][ix]+g->bz[iz+_C.idz][iy][ix]));
}
}
if(MASTER) fprintf(_FLOG,"..INFO(configUser): MHD/PROBLEM24 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/problem24B: Kelvin-Helmholtz instability
*/
>SIMULATION I/O ----------------------------------------------------------------------------------
NEW ./NIRLAST >mode:{NEW,RES,MOD,ANA},fname
1000000 1.59e+02 >mod_max,time_max
0100 02000 00500 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.00000000e+00 2.00000000e+01 0128 >lo[0],up[0],dim[0]
-1.00000000e+01 1.00000000e+01 0128 >lo[1],up[1],dim[1]
0.00000000e+00 1.00000000e+00 0000 >lo[2],up[2],dim[2]
SLOCK 1 1 0 >partitioning_type:{SFC,BLOCK # # #}
>BOUNDARY CONDITIONS -----------------------------------------------------------------------------
PPMMOO >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.00000e+01 0.00000e+00 0.00000e+00 >param[0-2] / Ma
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] /
0 0 0 0 0 0 0 0 >flag[0-7] /
>SOLVER SPECIFICATIONS ---------------------------------------------------------------------------
RK3 HLLD CCT 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 2.00e+00 1.00e+00 1.00e+00 >eos,gamma,poly_const,temperature
N >heatloss
0 >tracer
0 >fields
N N N 1.0 0.0 >species,reac,energy_reac,mmw,mean_ionization
>>END INPUT --------------------------------------------------------------------------------------
pro readCSV_MHD_24B,FILE,M,V
;
; NIRVANA: CAIVS: IDL
;
; Version: 4
;
; Reads CSV-like data files (#.csv files) produced by testproblem
; MHD/problem24B.
;
; Copyright (C) 2020 Udo Ziegler
;
; @param FILE full name of data file
; @param M metadata header
; @param V data array
;
; @author Udo Ziegler, AIP
;-------------------------------------------------------------------------------
;
openr,1,FILE
print
print,'IMPORTING FILE ',FILE,' ...'
; METADATA ---------------------------------------------------------------------
; BASIC ATTRIBUTES
str=''
ncols=0L
nrows=0L
readf,1,format='(A6,I3,A7,I9)',str,ncols,str,nrows
; COLUMN LABELS
label=''
readf,1,label
; ASSIGN OUTPUT HEADER ---------------------------------------------------------
M={ ncols:0L,nrows:0L,label:'' }
M.ncols=ncols
M.nrows=nrows
M.label=label
print
str=strcompress(string(M.ncols),remove_all=1)
print,' #cols: ',str
str=strcompress(string(M.nrows),remove_all=1)
print,'#rows: ',str
; READ DATA --------------------------------------------------------------------
VR=dblarr(ncols)
V=dblarr(nrows,ncols)
for ir=0,nrows-1 do begin
readf,1,VR
V(ir,*)=VR
endfor
print
print,'... IMPORT COMPLETED.'
;
close,1
;
return
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment