Commit 34bc6724 authored by m.tavelli's avatar m.tavelli

In debug, added tecplot files

parent 7034224d
......@@ -97,6 +97,7 @@ void FOCCZ4::FOCCZ4Solver_ADERDG::boundaryValues(const double* const x,const dou
std::memset(stateOut, 0, nVar * sizeof(double));
std::memset(fluxOut , 0, nVar * sizeof(double));
std::copy_n(stateIn,nVar,stateOut);
for(int dd=0; dd<nDim; dd++) F[dd] = Fs[dd];
for(int i=0; i < basisSize; i++) { // i == time
......@@ -109,7 +110,7 @@ void FOCCZ4::FOCCZ4Solver_ADERDG::boundaryValues(const double* const x,const dou
flux(Qgp, F);
for(int m=0; m < nVar; m++) {
stateOut[m] += weight * Qgp[m];
fluxOut[m] += weight * Fs[direction][m];
//fluxOut[m] += weight * Fs[direction][m];
}
}
}
......
......@@ -68,7 +68,7 @@ RECURSIVE SUBROUTINE InitialData(xGP, tGp, Q)
!
r = SQRT( xGP(1)**2 + xGP(2)**2 )
!
V0(60) = 0.001*EXP(-0.5*( (xGP(1)-2.0)**2+xGP(2)**2+0*xGP(3)**2)/1.0**2 ) + 0.001*EXP(-0.5*( (xGP(1)+2.0)**2+xGP(2)**2+0*xGP(3)**2)/1.0**2 )
!V0(60) = 0.001*EXP(-0.5*( (xGP(1)-2.0)**2+xGP(2)**2+0*xGP(3)**2)/1.0**2 ) + 0.001*EXP(-0.5*( (xGP(1)+2.0)**2+xGP(2)**2+0*xGP(3)**2)/1.0**2 )
IF( r>5.0 .AND. r<10. ) THEN
V0(61) = -0.2*xGP(2)*( 10.0 - r )/5.0
V0(62) = +0.2*xGP(1)*( 10.0 - r )/5.0
......@@ -79,7 +79,7 @@ RECURSIVE SUBROUTINE InitialData(xGP, tGp, Q)
V0(61:62) = 0.0
ENDIF
V0(63) = 0.0
V0(64) = 0.5*0.001*EXP(-0.5*( (xGP(1)-2.0)**2+xGP(2)**2+0*xGP(3)**2)/1.0**2 ) + 0.5*0.001*EXP(-0.5*( (xGP(1)+2.0)**2+xGP(2)**2+0*xGP(3)**2)/1.0**2 )
!V0(64) = 0.5*0.001*EXP(-0.5*( (xGP(1)-2.0)**2+xGP(2)**2+0*xGP(3)**2)/1.0**2 ) + 0.5*0.001*EXP(-0.5*( (xGP(1)+2.0)**2+xGP(2)**2+0*xGP(3)**2)/1.0**2 )
!
CASE DEFAULT
PRINT *, 'NOT IMPLEMENTED'
......
......@@ -36,14 +36,16 @@ RECURSIVE SUBROUTINE PDEFlux(f,g,hz,Q)
END SUBROUTINE PDEFlux
RECURSIVE SUBROUTINE PDENCP(BgradQ,Q,gradQ)
RECURSIVE SUBROUTINE PDENCP(BgradQ,Q,gradQin)
USE MainVariables, ONLY : nVar, nDim, EQN,nParam,d
IMPLICIT NONE
! 11. Oct 21:40: This was a matrix BGradQ(nVar, nDim) but is a vector in spaceTimePredictorNonlinear
REAL, INTENT(OUT) :: BgradQ(nVar)
REAL, INTENT(IN) :: gradQ(nVar, 3)
REAL, INTENT(IN) :: gradQin(nVar, 3)
REAL, INTENT(IN) :: Q(nVar)
REAL :: gradQ(nVar, 3)
! Argument list
REAL :: par(nParam)
! Local variables
......@@ -90,6 +92,7 @@ RECURSIVE SUBROUTINE PDENCP(BgradQ,Q,gradQ)
#endif
!
BgradQ = 0.0
!return
!
#ifdef EULER
!
......@@ -265,12 +268,15 @@ BgradQ(7) = -6.0*EQN%c0**2*Q(2)/Q(1)*gradQ(8,1)
!
#if defined(CCZ4EINSTEIN) || defined(CCZ4GRMHD) || defined(CCZ4GRHD) || defined(CCZ4GRGPR)
!
Qx = gradQ(:,1)
Qy = gradQ(:,2)
Qx = gradQin(:,1)
Qy = gradQin(:,2)
IF(nDim.eq.2) THEN
Qz = 0.0
gradQ(:,1:2)=gradQin(:,1:2)
gradQ(:,3) =0.0
ELSE
Qz = gradQ(:,3)
Qz = gradQin(:,3)
gradQ=gradQin
ENDIF
k1 = EQN%CCZ4k1
......@@ -2301,13 +2307,14 @@ END SUBROUTINE PDEMatrixB
!END SUBROUTINE getExactSolution
SUBROUTINE PDEFusedSrcNCP(Src_BgradQ,Q,gradQ)
USE MainVariables, ONLY : nVar, nParam, d, EQN
RECURSIVE SUBROUTINE PDEFusedSrcNCP(Src_BgradQ,Q,gradQin)
USE MainVariables, ONLY : nVar, nParam, d, EQN, nDim
IMPLICIT NONE
! Argument list
REAL, INTENT(IN) :: Q(nVar), gradQ(nVar,d)
REAL, INTENT(IN) :: Q(nVar), gradQin(nVar,d)
REAL, INTENT(OUT) :: Src_BgradQ(nVar)
! Local variables
REAL :: gradQ(nVar,d)
REAL :: par(nParam), time=0.0
INTEGER :: i,j,k,l,m,n,ip,iq,ii,jj,kk,ll,mm,iErr,count
REAL :: p, irho, lam, mu
......@@ -2352,7 +2359,7 @@ SUBROUTINE PDEFusedSrcNCP(Src_BgradQ,Q,gradQ)
! Matter contributions
REAL :: sm(3), Sij(3,3), Sij_contr(3,3), sm_contr(3), S, tau, dens, bv_cov(3), sqrdet
REAL :: srctraceK, srcTheta
REAL :: SrcK(3,3), SrcGhat(3)
REAL :: SrcK(3,3), SrcGhat(3),mytmp1(3,3,3,3), mytmp2(3,3,3,3)
REAL, PARAMETER :: Pi = ACOS(-1.0)
!
......@@ -2360,9 +2367,16 @@ SUBROUTINE PDEFusedSrcNCP(Src_BgradQ,Q,gradQ)
!
#if defined(CCZ4EINSTEIN) || defined(CCZ4GRHD) || defined(CCZ4GRMHD) || defined(CCZ4GRGPR)
!
Qx = gradQ(:,1)
Qy = gradQ(:,2)
Qz = gradQ(:,3)
Qx = gradQin(:,1)
Qy = gradQin(:,2)
IF(nDim.eq.2) THEN
Qz = 0.0
gradQ(:,1:2)=gradQin(:,1:2)
gradQ(:,3) =0.0
ELSE
Qz = gradQin(:,3)
gradQ=gradQin
ENDIF
k1 = EQN%CCZ4k1
k2 = EQN%CCZ4k2
......@@ -2577,13 +2591,25 @@ SUBROUTINE PDEFusedSrcNCP(Src_BgradQ,Q,gradQ)
Christoffel_kind1(i,j,k) = DD(k,i,j)+DD(j,i,k)-DD(i,j,k) ! this definition seems to work !
DO l = 1, 3
Christoffel_tilde(i,j,k) = Christoffel_tilde(i,j,k) + g_contr(k,l)*( DD(i,j,l)+DD(j,i,l)-DD(l,i,j) )
Christoffel(i,j,k) = Christoffel(i,j,k) + g_contr(k,l)*( DD(i,j,l)+DD(j,i,l)-DD(l,i,j) ) -g_contr(k,l)*( g_cov(j,l)*PP(i)+g_cov(i,l)*PP(j)-g_cov(i,j)*PP(l) )
mytmp1(i,j,k,l) = DD(i,j,l)+DD(j,i,l)-DD(l,i,j)
Christoffel(i,j,k) = Christoffel(i,j,k) + g_contr(k,l)*( mytmp1(i,j,k,l) )
mytmp2(i,j,k,l)=( g_cov(j,l)*PP(i)+g_cov(i,l)*PP(j)-g_cov(i,j)*PP(l) )
Christoffel(i,j,k) = Christoffel(i,j,k) -g_contr(k,l)*mytmp2(i,j,k,l)
!Gtilde(i) = Gtilde(i)+2*g_contr(i,j)*g_contr(k,l)*DD(l,j,k)
!PRINT *, Christoffel(i,j,k) ,I,J,K
ENDDO
ENDDO
ENDDO
ENDDO
!
!print *, 'g_contr',g_contr
!print *, 'g_cov',g_cov
!print *, 'DD',DD
!print *, 'PP',PP
print *, 'mytmp1',mytmp1
print *, 'mytmp2',mytmp2
print *, 'Christoffel',Christoffel
pause
DO l = 1, 3
DO j = 1, 3
DO i = 1, 3
......@@ -2737,6 +2763,23 @@ SUBROUTINE PDEFusedSrcNCP(Src_BgradQ,Q,gradQ)
!
dtTraceK = - nablanablaalphaNCP - nablanablaalphaSrc + alpha*( RPlusTwoNablaZNCP + RPlusTwoNablaZSrc + traceK**2 - 2*Theta*traceK ) - 3*alpha*k1*(1+k2)*Theta + SUM(beta(:)*dtraceK(:))
!
!print *, 'dtTraceK=', dtTraceK
!print *, 'nablanablaalphaNCP=', nablanablaalphaNCP
!print *, 'nablanablaalphaSrc=', nablanablaalphaSrc
!print *, 'RPlusTwoNablaZNCP=', RPlusTwoNablaZNCP
!print *, 'RPlusTwoNablaZSrc=', RPlusTwoNablaZSrc
!print *, 'traceK=', traceK
!print *, 'dtraceK=', dtraceK
!print *, 'Theta=', Theta
!print *, 'Christoffel',Christoffel
print *, 'Christoffel',Christoffel
!print *, 'Gtilde',Gtilde
!print *, 'Christoffel_tilde',Christoffel_tilde
!print *, 'dChristoffelsrc',dChristoffelsrc
!print *, 'dChristoffel_tildesrc',dChristoffel_tildesrc
!print *, 'Z',Z
pause
traceB = BB(1,1) + BB(2,2) + BB(3,3)
dtphi = beta(1)*PP(1) + beta(2)*PP(2) + beta(3)*PP(3) + 1./3.*alpha*traceK - 1./3.*traceB
dtalpha = -alpha*fa*(traceK-K0-c*2*Theta) + beta(1)*AA(1) + beta(2)*AA(2) + beta(3)*AA(3)
......
This diff is collapsed.
RECURSIVE SUBROUTINE ElementCallTECPLOTPLOTTER(wh_in,lx0,ldx,limiter)
USE TECPLOTPLOTTERmod
USE MainVariables, only : nVar, nDim
implicit none
REAL, INTENT(IN) :: wh_in(nVar,nDOFm),lx0(nDim),ldx(nDim)
real :: wh(nVar,nDOFm)
real :: lx0_3(3),ldx_3(3)
integer :: limiter
wh=wh_in
lx0_3=0.
ldx_3=0.
lx0_3(1:nDim)=lx0
ldx_3(1:nDim)=ldx
CALL ElementTECPLOTPLOTTER(wh,lx0_3,ldx_3,limiter)
END SUBROUTINE ElementCallTECPLOTPLOTTER
RECURSIVE SUBROUTINE ElementCallTECPLOTFVPLOTTER(wh,lx0,ldx,limiter)
USE TECPLOTPLOTTERmod
USE MainVariables, only : nVar, nDim
implicit none
REAL, INTENT(IN) :: wh(nVar,(nSubLim+2*nSubLim_GL)**nDIM),lx0(nDim),ldx(nDim)
real :: lx0_3(3),ldx_3(3)
integer :: limiter,i,j,k,Stencil
lx0_3(1:nDim)=lx0
ldx_3(1:nDim)=ldx
CALL ElementTECPLOTPLOTTER_FV(wh,lx0_3,ldx_3,limiter)
END SUBROUTINE ElementCallTECPLOTFVPLOTTER
RECURSIVE SUBROUTINE InitializeTECPLOTPLOTTER(time)
USE TECPLOTPLOTTERmod
implicit none
real :: time
CALL InitTECPLOTPLOTTER(time)
END SUBROUTINE InitializeTECPLOTPLOTTER
RECURSIVE SUBROUTINE FinishTECPLOTPLOTTER(Myrank)
USE TECPLOTPLOTTERmod
implicit none
integer :: Myrank
CALL FinalizeTECPLOTPLOTTER(Myrank)
END SUBROUTINE FinishTECPLOTPLOTTER
#ifndef __EXAHYPE_USER_TECINT__
#define __EXAHYPE_USER_TECINT__
// Fortran functions:
extern "C" {
//void inittecplot_(const int* N_in, const int* M_in, const int* basisSize, const int* Ghostlayers);
void elementcalltecplotplotter_(const double *wh, const double* lx0, const double* ldx, const int* limiter);
void elementcalltecplotaderdgplotter_(const double *wh, const double* lx0, const double* ldx, const int* limiter);
void elementcalltecplotfvplotter_(const double *wh, const double* lx0, const double* ldx, const int* limiter);
void finishtecplotplotter_(const int* Myrank);
void initializetecplotplotter_(const double* time);
}/* extern "C" */
#endif /* __EXAHYPE_USER_PDE__ */
// This file was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
//
// ========================
// www.exahype.eu
// ========================
#include "TecplotWriter.h"
#include "TECPLOTinterface.h"
#include "tarch/parallel/Node.h"
#include "tarch/parallel/NodePool.h"
#include <stdio.h>
#include "FOCCZ4Solver_FV.h"
#include "FOCCZ4Solver_ADERDG.h"
#include "kernels/GaussLegendreBasis.h"
#include "kernels/KernelUtils.h"
#include "peano/utils/Loop.h"
#include "tarch/la/VectorOperations.h"
#include <algorithm>
#include <iomanip>
FOCCZ4::TecplotWriter::TecplotWriter() : exahype::plotters::LimitingADERDG2UserDefined::LimitingADERDG2UserDefined(){
// @TODO Please insert your code here.
}
void FOCCZ4::TecplotWriter::plotADERDGPatch(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch, double* const u,
double timeStamp) {
// @TODO Please insert your code here.
plotForADERSolver = 0;
elementcalltecplotplotter_(u,&offsetOfPatch[0],&sizeOfPatch[0],&plotForADERSolver);
}
void FOCCZ4::TecplotWriter::plotFiniteVolumesPatch(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch, double* const u,
double timeStamp) {
// @TODO Please insert your code here.
plotForADERSolver = 1;
elementcalltecplotfvplotter_(u, &offsetOfPatch[0], &sizeOfPatch[0], &plotForADERSolver);
}
void FOCCZ4::TecplotWriter::startPlotting( double time) {
// @TODO Please insert your code here.
mpirank = tarch::parallel::Node::getInstance().getRank();
initializetecplotplotter_(&time);
}
void FOCCZ4::TecplotWriter::finishPlotting() {
// @TODO Please insert your code here.
finishtecplotplotter_(&mpirank);
}
\ No newline at end of file
// This file was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
// ========================
// www.exahype.eu
// ========================
#ifndef TecplotWriter_CLASS_HEADER_
#define TecplotWriter_CLASS_HEADER_
#include "exahype/plotters/LimitingADERDG2UserDefined.h"
#include "exahype/plotters/ADERDG2UserDefined.h"
#include "exahype/plotters/ascii/MultipleReductionsWriter.h"
#include "exahype/solvers/LimitingADERDGSolver.h"
#include "FOCCZ4Solver_ADERDG_Variables.h"
#include "FOCCZ4Solver_ADERDG.h"
#include "FOCCZ4Solver_FV.h"
namespace FOCCZ4 {
class TecplotWriter;
}
class FOCCZ4::TecplotWriter : public exahype::plotters::LimitingADERDG2UserDefined {
private:
int plotForADERSolver;
int plotForFVSolver;
int mpirank;
public:
/**
* Constructor.
*
* \note ExaHyPE does not increment file counters for
* you if you use user defined plotting. You have
* to declare and manage such member variables yourself.
*/
TecplotWriter();
/**
* This method is invoked every time an ADER-DG cell
* is touched by the plotting plotter.
*
* \note Use the protected variables _order, _variables to
* determine the size of u.
* The array u has the size _variables * (_order+1)^DIMENSIONS.
*
* \param[in] offsetOfPatch the offset of the cell/patch.
* \param[in] sizeOfPatch the offset of the cell/patch.
* \param[in] u the degrees of freedom "living" inside of the patch.
*/
void plotADERDGPatch(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch, double* const u,
double timeStamp) override;
/**
* This method is invoked every time a Finite Volume cell
* is touched by the plotting plotter.
*
* \note Use the protected variables _order, _variables to
* determine the size of u.
* The array u has the size _variables * (_order+1)^DIMENSIONS.
*
* \param[in] offsetOfPatch the offset of the cell/patch.
* \param[in] sizeOfPatch the offset of the cell/patch.
* \param[in] u the degrees of freedom "living" inside of the patch.
*/
void plotFiniteVolumesPatch(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch, double* const u,
double timeStamp) override;
/**
* This method is called at the beginning of the plotting.
* You can use it to reset member variables, e.g., those
* used for calculations, or to increment file counters.
*
* \param[in] time a characteristic solver time stamp.
* Usually the global minimum.
*/
void startPlotting( double time) override;
/**
* This method is called at the end of the plotting.
* You can use it to reset member variables, finalise calculations (compute square roots etc.),
* or to increment file counters
*/
void finishPlotting() override;
};
#endif /* TecplotWriter_CLASS_HEADER_ */
\ No newline at end of file
! Tools.f90
RECURSIVE SUBROUTINE InitTECPLOT(N_in,M_in,SubLim_in,Ghostlayers_in)
USE TECPLOTPLOTTERmod
implicit none
INTEGER :: N_in,M_in,SubLim_in,Ghostlayers_in
CALL SetMainParameters(N_in,M_in,SubLim_in,Ghostlayers_in)
END SUBROUTINE InitTECPLOT
!RECURSIVE SUBROUTINE InitTECPLOT(N_in,SubLim_in,Ghostlayers_in)
! USE, INTRINSIC :: ISO_C_BINDING
! USE TECPLOTPLOTTERmod !, only : SetMainParameters
! implicit none
! INTEGER :: N_in,SubLim_in,Ghostlayers_in
! CALL SetMainParameters(N_in,SubLim_in,Ghostlayers_in)
!END SUBROUTINE InitTECPLOT
RECURSIVE SUBROUTINE getNumericalSolution(V,Q)
USE MainVariables, ONLY: nVar
IMPLICIT NONE
REAL :: V(nVar), Q(nVar)
INTEGER :: iErr
!
CALL PDECons2Prim(V,Q,iErr)
!
END SUBROUTINE getNumericalSolution
RECURSIVE SUBROUTINE getExactSolution(x,timeStamp,V)
USE MainVariables, ONLY: nAux,nVar,nDim
IMPLICIT NONE
REAL :: V(nVar), Q(nVar), x(nDim), timeStamp
INTEGER :: iErr
!
call InitialData(x, timeStamp, Q)
CALL PDECons2Prim(V,Q,iErr)
!
END SUBROUTINE getExactSolution
// Tools.h
// Fortran functions:
extern "C" {
//
void getnumericalsolution_(double* V,double* Q);
void getexactsolution_(double* x,double* timestep,double* V);
void inittecplot_(const int* N_in,const int* M_in, const int* basisSize, const int* Ghostlayers);
//void inittecplot_(const int* N_in,const int* M_in);
//
}/* extern "C" */
\ No newline at end of file
{
"project_name": "FOCCZ4",
"compiler_flags": "-DCCZ4EINSTEIN -DGLMROT",
"paths": {
"peano_kernel_path": "./Peano",
"exahype_path": "./ExaHyPE",
"output_directory": "./ApplicationExamples/FOCCZ4/FOCCZ4",
"log_file": "whatever.log"
},
"architecture": "skx",
"computational_domain": {
"dimension": 3,
"end_time": 3.0,
"offset": [
-20.0,
-20.0,
-20.0
],
"width": [
40.0,
40.0,
40.0
]
},
"shared_memory": {
"cores": 1,
"properties_file": "sharedmemory.properties",
"autotuning_strategy": "dummy",
"background_job_consumers": 1
},
"distributed_memory": {
"timeout": 6000,
"load_balancing_type": "static",
"buffer_size": 6400,
"load_balancing_strategy": "hotspot",
"node_pool_strategy": "fair",
"ranks_per_node": 10
},
"optimisation": {
"fuse_algorithmic_steps": "all",
"fuse_algorithmic_steps_rerun_factor": 0.99,
"fuse_algorithmic_steps_diffusion_factor": 0.99,
"spawn_predictor_as_background_thread": true,
"spawn_update_as_background_thread": true,
"spawn_amr_background_threads": true,
"disable_vertex_exchange_in_time_steps": true,
"time_step_batch_factor": 0.0,
"disable_metadata_exchange_in_batched_time_steps": false,
"double_compression": 0.0,
"spawn_double_compression_as_background_thread": true
},
"solvers": [
{
"type": "Limiting-ADER-DG",
"name": "FOCCZ4Solver",
"order": 3,
"maximum_mesh_size": 10.0,
"maximum_mesh_depth": 0,
"time_stepping": "global",
"aderdg_kernel": {
"language": "C",
"nonlinear": true,
"terms": [
"flux",
"ncp",
"source"
],
"space_time_predictor": {},
"optimised_terms": [
"fusedsource"
],
"optimised_kernel_debugging": [],
"implementation": "optimised"
},
"point_sources": 0,
"limiter": {
"dmp_observables": 3,
"dmp_relaxation_parameter": 1e+3,
"dmp_difference_scaling": 1e+4,
"patch_size": "max",
"implementation": "generic"
},
"fv_kernel": {
"language": "C",
"terms": [
"flux",
"ncp",
"source"
],
"scheme": "musclhancock",
"slope_limiter" : "minmod",
"implementation": "generic"
},
"variables": [
{
"name": "MYQ",
"multiplicity": 96
}
],
"parameters": {
"reference": "CCZ4MinkowskiSrc"
},
"plotters": [
{
"type": "user::defined",
"name": "TecplotWriter",
"time": 0.0,
"repeat": 0.05,
"output": "./output/tecplot",
"variables": 96
}
]
}
]
}
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