Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature/hevi no mass matrix lumping #381

Merged
merged 2 commits into from
Dec 28, 2024
Merged

Feature/hevi no mass matrix lumping #381

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0556081
Modify Makefile
ywkawai Dec 28, 2024
8c5f93f
Modify band matrix for HEVI to support a case without mass matrix lum…
ywkawai Dec 28, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions FElib/src/Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -286,9 +286,9 @@ info:

prep_scalelib:
echo "Prepair libscale_sub.a and include files for FE-Project .."
$(MAKE) $(LIBDIR)/libscale_sub.a
$(MAKE) $(LIBDIR)/libscale_sub$(POSTFIX).a

$(LIBDIR)/libscale_sub.a: $(SCALE)/lib/libscale.a
$(LIBDIR)/libscale_sub$(POSTFIX).a: $(SCALE)/lib/$(SCALE_LIBNAME)
mkdir -p $(LIBDIR)
$(AR) -x $< $(OBJS_NAME_SCALELIB)
$(AR) $(ARFLAGS) $(@F) $(OBJS_NAME_SCALELIB)
Expand Down
7 changes: 3 additions & 4 deletions FElib/src/fluid_dyn_solver/scale_atm_dyn_dgm_nonhydro3d_rhot_hevi.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -283,12 +283,11 @@ subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_cal_vi( &
lmesh, elem, lmesh2D, elem2D ) ! (in)

use scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common, only: &
vi_gen_vmap => atm_dyn_dgm_nonhydro3d_rhot_hevi_common_gen_vmap, &
vi_eval_Ax => atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_Ax, &
vi_gen_vmap => atm_dyn_dgm_nonhydro3d_rhot_hevi_common_gen_vmap, &
vi_eval_Ax => atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_Ax, &
vi_eval_Ax_uv => atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_Ax_uv, &
vi_construct_matbnd => atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd, &
vi_construct_matbnd => atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd, &
vi_construct_matbnd_uv => atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv

implicit none

class(LocalMesh3D), intent(in) :: lmesh
Expand Down
67 changes: 34 additions & 33 deletions FElib/src/fluid_dyn_solver/scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -502,7 +502,7 @@ subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd( &
integer :: Colmask(elem%Nnode_v)
real(RP) :: Id(elem%Nnode_v,elem%Nnode_v)
real(RP) :: Dd(elem%Nnode_v)
real(RP) :: tmp1
real(RP) :: tmp1(elem%Nnode_v)
real(RP) :: fac
integer :: ij, v1, v2, pv1, pv2, g_kj, g_kjp1, g_kjm1, pb1
logical :: bc_flag
Expand Down Expand Up @@ -628,64 +628,65 @@ subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd( &
fp2 = elem%Nfp_h * elem%Nfaces_h + (f2-1)*elem%Nfp_v + ij

!--
tmp1 = fac * elem%Lift(FmV,fp) * lmesh%Fscale(fp,ke) &
tmp1(:) = fac * elem%Lift(Colmask(:),fp) * lmesh%Fscale(fp,ke) &
* max( alph(fp,ke_z), alph(fp2,ke_z2) )
if (bc_flag) then
PmatD(pv1,pv1,MOMZ_VID_LC,MOMZ_VID_LC) = PmatD(pv1,pv1,MOMZ_VID_LC,MOMZ_VID_LC) + 2.0_RP * tmp1
PmatD(:,pv1,MOMZ_VID_LC,MOMZ_VID_LC) = PmatD(:,pv1,MOMZ_VID_LC,MOMZ_VID_LC) + 2.0_RP * tmp1(:)
else
do v=1, 3
PmatD(pv1,pv1,v,v) = PmatD(pv1,pv1,v,v) + tmp1
PmatD(:,pv1,v,v) = PmatD(:,pv1,v,v) + tmp1(:)
if (f1 == 1) then
PmatL(pv1,pv2,v,v) = - tmp1
PmatL(:,pv2,v,v) = - tmp1(:)
else
PmatU(pv1,pv2,v,v) = - tmp1
PmatU(:,pv2,v,v) = - tmp1(:)
end if
end do
end if

!--
tmp1 = fac * elem%Lift(FmV,fp) * lmesh%Fscale(fp,ke) * nz(fp,ke_z)
tmp1(:) = fac * elem%Lift(Colmask(:),fp) * lmesh%Fscale(fp,ke) * nz(fp,ke_z)

if (bc_flag) then
PmatD(pv1,pv1,DENS_VID_LC,MOMZ_VID_LC) = PmatD(pv1,pv1,DENS_VID_LC,MOMZ_VID_LC) - 2.0_RP * tmp1
PmatD(pv1,pv1,RHOT_VID_LC,MOMZ_VID_LC) = PmatD(pv1,pv1,RHOT_VID_LC,MOMZ_VID_LC) - 2.0_RP * tmp1 * POT0(pv1,ke_z,ij)
PmatD(pv1,pv1,RHOT_VID_LC,DENS_VID_LC) = PmatD(pv1,pv1,RHOT_VID_LC,DENS_VID_LC) + 2.0_RP * tmp1 * POT0(pv1,ke_z,ij) * WT0(pv1,ke_z,ij)
PmatD(pv1,pv1,RHOT_VID_LC,RHOT_VID_LC) = PmatD(pv1,pv1,RHOT_VID_LC,RHOT_VID_LC) - 2.0_RP * tmp1 * WT0(pv1,ke_z,ij)
PmatD(:,pv1,DENS_VID_LC,MOMZ_VID_LC) = PmatD(:,pv1,DENS_VID_LC,MOMZ_VID_LC) - 2.0_RP * tmp1(:)
PmatD(:,pv1,RHOT_VID_LC,MOMZ_VID_LC) = PmatD(:,pv1,RHOT_VID_LC,MOMZ_VID_LC) - 2.0_RP * tmp1(:) * POT0(pv1,ke_z,ij)
PmatD(:,pv1,RHOT_VID_LC,DENS_VID_LC) = PmatD(:,pv1,RHOT_VID_LC,DENS_VID_LC) + 2.0_RP * tmp1(:) * POT0(pv1,ke_z,ij) * WT0(pv1,ke_z,ij)
PmatD(:,pv1,RHOT_VID_LC,RHOT_VID_LC) = PmatD(:,pv1,RHOT_VID_LC,RHOT_VID_LC) - 2.0_RP * tmp1(:) * WT0(pv1,ke_z,ij)
else
PmatD(pv1,pv1,DENS_VID_LC,MOMZ_VID_LC) = PmatD(pv1,pv1,DENS_VID_LC,MOMZ_VID_LC) - tmp1
PmatD(:,pv1,DENS_VID_LC,MOMZ_VID_LC) = PmatD(:,pv1,DENS_VID_LC,MOMZ_VID_LC) - tmp1(:)

! PmatD(pv1,pv1,MOMZ_VID_LC,MOMZ_VID_LC) = PmatD(pv1,pv1,MOMZ_VID_LC,MOMZ_VID_LC) - tmp1 * 2.0_RP * W0(pv1,ke_z,ij) [ <- d_MOMZ ( MOMZ x MOMZ / DENS ) ]
! PmatD(pv1,pv1,MOMZ_VID_LC,DENS_VID_LC) = PmatD(pv1,pv1,MOMZ_VID_LC,DENS_VID_LC) + tmp1 * W0(pv1,ke_z,ij)**2 [ <- d_DENS ( MOMZ x MOMZ / DENS ) ]
PmatD(pv1,pv1,MOMZ_VID_LC,RHOT_VID_LC) = PmatD(pv1,pv1,MOMZ_VID_LC,RHOT_VID_LC) - tmp1 * DPDRHOT0(pv1,ke_z,ij)
PmatD(:,pv1,MOMZ_VID_LC,RHOT_VID_LC) = PmatD(:,pv1,MOMZ_VID_LC,RHOT_VID_LC) - tmp1(:) * DPDRHOT0(pv1,ke_z,ij)

PmatD(pv1,pv1,RHOT_VID_LC,MOMZ_VID_LC) = PmatD(pv1,pv1,RHOT_VID_LC,MOMZ_VID_LC) - tmp1 * POT0(pv1,ke_z,ij)
PmatD(pv1,pv1,RHOT_VID_LC,DENS_VID_LC) = PmatD(pv1,pv1,RHOT_VID_LC,DENS_VID_LC) + tmp1 * POT0(pv1,ke_z,ij) * WT0(pv1,ke_z,ij)
PmatD(pv1,pv1,RHOT_VID_LC,RHOT_VID_LC) = PmatD(pv1,pv1,RHOT_VID_LC,RHOT_VID_LC) - tmp1 * WT0(pv1,ke_z,ij)
PmatD(:,pv1,RHOT_VID_LC,MOMZ_VID_LC) = PmatD(:,pv1,RHOT_VID_LC,MOMZ_VID_LC) - tmp1(:) * POT0(pv1,ke_z,ij)
PmatD(:,pv1,RHOT_VID_LC,DENS_VID_LC) = PmatD(:,pv1,RHOT_VID_LC,DENS_VID_LC) + tmp1(:) * POT0(pv1,ke_z,ij) * WT0(pv1,ke_z,ij)
PmatD(:,pv1,RHOT_VID_LC,RHOT_VID_LC) = PmatD(:,pv1,RHOT_VID_LC,RHOT_VID_LC) - tmp1(:) * WT0(pv1,ke_z,ij)

if (f1 == 1) then
PmatL(pv1,pv2,DENS_VID_LC,MOMZ_VID_LC) = + tmp1
! PmatL(pv1,pv2,DENS_VID_LC,MOMZ_VID_LC) = + tmp1
PmatL(:,pv2,DENS_VID_LC,MOMZ_VID_LC) = + tmp1(:)

! PmatL(pv1,pv2,MOMZ_VID_LC,MOMZ_VID_LC) = PmatL(pv1,pv2,MOMZ_VID_LC,MOMZ_VID_LC) & [ <- d_MOMZ ( MOMZ x MOMZ / DENS ) ]
! + tmp1 * 2.0_RP * W0(pv2,ke_z2,ij)
! PmatL(pv1,pv2,MOMZ_VID_LC,DENS_VID_LC) = - tmp1 * W0(pv2,ke_z2,ij) * W0(pv2,ke_z2,ij) [ <- d_DENS ( MOMZ x MOMZ / DENS ) ]
PmatL(pv1,pv2,MOMZ_VID_LC,RHOT_VID_LC) = + tmp1 * DPDRHOT0(pv2,ke_z2,ij)
PmatL(:,pv2,MOMZ_VID_LC,RHOT_VID_LC) = + tmp1(:) * DPDRHOT0(pv2,ke_z2,ij)

PmatL(pv1,pv2,RHOT_VID_LC,MOMZ_VID_LC) = + tmp1 * POT0(pv2,ke_z2,ij)
PmatL(pv1,pv2,RHOT_VID_LC,DENS_VID_LC) = - tmp1 * POT0(pv2,ke_z2,ij) * WT0(pv2,ke_z2,ij)
PmatL(pv1,pv2,RHOT_VID_LC,RHOT_VID_LC) = PmatL(pv1,pv2,RHOT_VID_LC,RHOT_VID_LC) &
+ tmp1 * WT0(pv2,ke_z2,ij)
PmatL(:,pv2,RHOT_VID_LC,MOMZ_VID_LC) = + tmp1(:) * POT0(pv2,ke_z2,ij)
PmatL(:,pv2,RHOT_VID_LC,DENS_VID_LC) = - tmp1(:) * POT0(pv2,ke_z2,ij) * WT0(pv2,ke_z2,ij)
PmatL(:,pv2,RHOT_VID_LC,RHOT_VID_LC) = PmatL(:,pv2,RHOT_VID_LC,RHOT_VID_LC) &
+ tmp1(:) * WT0(pv2,ke_z2,ij)
else
PmatU(pv1,pv2,DENS_VID_LC,MOMZ_VID_LC) = + tmp1
PmatU(:,pv2,DENS_VID_LC,MOMZ_VID_LC) = + tmp1(:)

! PmatU(pv1,pv2,MOMZ_VID_LC,MOMZ_VID_LC) = PmatU(pv1,pv2,MOMZ_VID_LC,MOMZ_VID_LC ) & [ <- d_MOMZ ( MOMZ x MOMZ / DENS ) ]
! + tmp1 * 2.0_RP * W0(pv2,ke_z2,ij)
! PmatU(pv1,pv2,MOMZ_VID_LC,DENS_VID_LC) = - tmp1 * W0(pv2,ke_z2,ij) * W0(pv2,ke_z2,ij) [ <- d_DENS ( MOMZ x MOMZ / DENS ) ]
PmatU(pv1,pv2,MOMZ_VID_LC,RHOT_VID_LC) = + tmp1 * DPDRHOT0(pv2,ke_z2,ij)
PmatU(:,pv2,MOMZ_VID_LC,RHOT_VID_LC) = + tmp1(:) * DPDRHOT0(pv2,ke_z2,ij)

PmatU(pv1,pv2,RHOT_VID_LC,MOMZ_VID_LC) = + tmp1 * POT0(pv2,ke_z2,ij)
PmatU(pv1,pv2,RHOT_VID_LC,DENS_VID_LC) = - tmp1 * POT0(pv2,ke_z2,ij) * WT0(pv2,ke_z2,ij)
PmatU(pv1,pv2,RHOT_VID_LC,RHOT_VID_LC) = PmatU(pv1,pv2,RHOT_VID_LC,RHOT_VID_LC) &
+ tmp1 * WT0(pv2,ke_z2,ij)
PmatU(:,pv2,RHOT_VID_LC,MOMZ_VID_LC) = + tmp1(:) * POT0(pv2,ke_z2,ij)
PmatU(:,pv2,RHOT_VID_LC,DENS_VID_LC) = - tmp1(:) * POT0(pv2,ke_z2,ij) * WT0(pv2,ke_z2,ij)
PmatU(:,pv2,RHOT_VID_LC,RHOT_VID_LC) = PmatU(:,pv2,RHOT_VID_LC,RHOT_VID_LC) &
+ tmp1(:) * WT0(pv2,ke_z2,ij)
end if
end if
end do
Expand Down Expand Up @@ -777,7 +778,7 @@ subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv( &
integer :: Colmask(elem%Nnode_v)
real(RP) :: Id(elem%Nnode_v,elem%Nnode_v)
real(RP) :: Dd(elem%Nnode_v)
real(RP) :: tmp1
real(RP) :: tmp1(elem%Nnode_v)
real(RP) :: fac
integer :: ij, v1, v2, pv1, pv2, g_kj, g_kjp1, g_kjm1, pb1
logical :: bc_flag
Expand Down Expand Up @@ -846,15 +847,15 @@ subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv( &
fp2 = elem%Nfp_h * elem%Nfaces_h + (f2-1)*elem%Nfp_v + ij

!--
tmp1 = fac * elem%Lift(FmV,fp) * lmesh%Fscale(fp,ke) &
tmp1(:) = fac * elem%Lift(Colmask(:),fp) * lmesh%Fscale(fp,ke) &
* max( alph(fp,ke_z), alph(fp2,ke_z2) )
if (bc_flag) then
else
PmatD_uv(pv1,pv1) = PmatD_uv(pv1,pv1) + tmp1
PmatD_uv(:,pv1) = PmatD_uv(:,pv1) + tmp1(:)
if (f1 == 1) then
PmatL_uv(pv1,pv2) = - tmp1
PmatL_uv(:,pv2) = - tmp1(:)
else
PmatU_uv(pv1,pv2) = - tmp1
PmatU_uv(:,pv2) = - tmp1(:)
end if
end if
end do
Expand Down
58 changes: 31 additions & 27 deletions Mkinclude
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,45 +1,29 @@
.DEFAULT_GOAL := all

OBJ_DIR = ./.libs
INST_DIR = $(TOPDIR)

ifeq ($(origin BINDIR),undefined)
BINDIR = $(INST_DIR)/bin
endif

ifeq ($(origin LIBDIR),undefined)
LIBDIR = $(INST_DIR)/lib
endif

ifeq ($(origin MODDIR),undefined)
MODDIR = $(INST_DIR)/include
endif

DCUTILSDIR = $(TOPDIR)/dc_utils
DCUTILSDIR = $(TOPDIR)/dc_utils
SCALEFELIBDIR = $(TOPDIR)/FElib
GLOBALSWDIR = $(TOPDIR)/model/global_shallow_water
GLOBALSWDIR = $(TOPDIR)/model/global_shallow_water
ATMNONHYDRO2DDIR = $(TOPDIR)/model/atm_nonhydro2d
ATMNONHYDRO3DDIR = $(TOPDIR)/model/atm_nonhydro3d

SCALE_INCLUDE = -I$(MODDIR)/scalelib -I$(SCALE)/scalelib/include
SCALE_LIBS = -L$(LIBDIR) -lscale_sub -L$(SCALE)/lib -ldcutils

CONTRIB_LIBS = $(SCALE_LIBS)
CONTRIB_INCLUDE = $(SCALE_INCLUDE)

ifeq ($(SCALE_DEBUG),T)
FFLAGS = $(FFLAGS_DEBUG) $(FFLAGS_SYSDEPEND) -DDEBUG $(SCALE_INCLUDE)
FFLAGS = $(FFLAGS_DEBUG) $(FFLAGS_SYSDEPEND) -DDEBUG
CFLAGS = $(CFLAGS_DEBUG) $(CFLAGS_SYSDEPEND)
POSTFIX = _debug
else
FFLAGS = $(FFLAGS_FAST) $(FFLAGS_SYSDEPEND) $(SCALE_INCLUDE)
FFLAGS = $(FFLAGS_FAST) $(FFLAGS_SYSDEPEND)
CFLAGS = $(CFLAGS_FAST) $(CFLAGS_SYSDEPEND)
ifeq ($(SCALE_USE_AGGRESSIVEOPT),T)
FFLAGS += $(FFLAGS_AGGRESSIVE)
CFLAGS += $(CFLAGS_AGGRESSIVE)
endif
ifeq ($(SCALE_QUICKDEBUG),T)
FFLAGS += $(FFLAGS_QUICKDEBUG) -DQUICKDEBUG
POSTFIX = _quickdebug
else
POSTFIX =
endif
endif

Expand All @@ -51,10 +35,6 @@ ifeq ($(SCALE_USE_SINGLEFP),T)
FFLAGS += -DSINGLE
endif

ifeq ($(SCALE_USE_MASSCHECK),T)
FFLAGS += -DCHECK_MASS
endif

ifeq ($(SCALE_ENABLE_OPENMP),T)
FFLAGS += $(FFLAGS_OPENMP)
endif
Expand All @@ -67,6 +47,30 @@ ifeq ($(SCALE_DEVELOP),T)
FFLAGS += -DSCALE_DEVELOP
endif

# ----
ifeq ($(origin BINDIR),undefined)
BINDIR = $(INST_DIR)/bin
endif

ifeq ($(origin LIBDIR),undefined)
LIBDIR = $(INST_DIR)/lib
endif

ifeq ($(origin MODDIR),undefined)
MODDIR = $(INST_DIR)/include$(POSTFIX)
endif

# SCALE library
ifeq ($(origin SCALE_LIBNAME),undefined)
SCALE_LIBNAME=libscale.a
endif
SCALE_INCLUDE = -I$(MODDIR)/scalelib -I$(SCALE)/scalelib/include
SCALE_LIBS = -L$(LIBDIR) -lscale_sub$(POSTFIX) -L$(SCALE)/lib -ldcutils
FFLAGS += $(SCALE_INCLUDE)

CONTRIB_LIBS = $(SCALE_LIBS)
CONTRIB_INCLUDE = $(SCALE_INCLUDE)


# NetCDF library setting
# library location are inquired in order of: environment variable->Makedef.XXX->here
Expand Down
Loading